Line data Source code
1 : //# ClarkCleanLatModel.h: this defines ClarkCleanLatModel
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_CLARKCLEANLATMODEL_H
30 : #define SYNTHESIS_CLARKCLEANLATMODEL_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <synthesis/MeasurementEquations/LinearModel.h>
34 : #include <synthesis/MeasurementEquations/LinearEquation.h>
35 : #include <casacore/lattices/Lattices/Lattice.h>
36 : #include <casacore/casa/Arrays/IPosition.h>
37 : #include <synthesis/MeasurementEquations/Iterate.h>
38 : #include <casacore/casa/Logging/LogIO.h>
39 : #include <casacore/casa/Arrays/ArrayFwd.h>
40 :
41 : namespace casa { //# NAMESPACE CASA - BEGIN
42 :
43 : class ClarkCleanProgress;
44 : class LatConvEquation;
45 : class CCList;
46 :
47 : // <summary>
48 : // A Class for performing the Clark Clean Algorithm on Arrays
49 : // </summary>
50 :
51 : // <use visibility=export>
52 :
53 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
54 : // </reviewed>
55 :
56 : // <prerequisite>
57 : // <li> ResidualEquation/LatConvEquation
58 : // <li> LinearModel/LinearEquation Paradigm
59 : // </prerequisite>
60 : //
61 : // <etymology>
62 : // This class is called ClarkCleanLatModel because thats the algorithm it uses
63 : // deconvolve the lattice-based model.
64 : // </etymology>
65 : //
66 : // <synopsis>
67 : // This class is used to perform the Clark Clean Algorithm on an
68 : // Array. Only the deconvolved model of the sky are directly stored by this
69 : // class. The point spread function (psf) and convolved (dirty) image are
70 : // stored in a companion class which is must be derived from
71 : // ResidualEquation.
72 : //
73 : // The cleaning works like this. The user constructs a ClarkCleanLatModel by
74 : // specifying an initial model of the sky. This can by be
75 : // one,two,three... dimensional depending on the dimension of the psf (see
76 : // below). The user then constructs a class which implements the forward
77 : // equation between the model and the dirty image. Typically this will be
78 : // the ConvolutionEquation class, although any class which has a
79 : // ResidualEquation interface will be work (but perhaps very slowly, as the
80 : // ConvolutionEquation class has member functions optimised for cleaning)
81 : //
82 : // The user then calls the solve() function (with the appropriate equation
83 : // class as an arguement), and this class will perform the Clark clean.
84 : // The various clean parameters are set (prior to calling solve) using the
85 : // functions derived from the Iterate class, in particular setGain(),
86 : // setNumberIterations() & setThreshold() (to set a flux limit).
87 : //
88 : // The solve() function does not return either the deconvolved model or the
89 : // residuals. The solved model can be obtained using the getModel() function
90 : // (derived from ArrayModel()) and the residual can be obtained using the
91 : // residual() member function of the Convolution/Residual Equation Class.
92 : //
93 : // The size and shape of the model used in this class MUST be the same as
94 : // the convolved data (Dirty Image), stored in the companion
95 : // ResidualEquation Class. However the model (and convolved data) can have
96 : // more dimensions than the psf, as well as a different size (either larger
97 : // or smaller). When the dimensionality is different the cleaning is done
98 : // independendtly in each "plane" of the model. (Note this has not
99 : // been implemented yet but is relatively simple to do if necessary).
100 : //
101 : // This multi-dimensionalty is exploited when cleaning arrays of
102 : // StokesVectors. Here the casacore::Array of StokesVectors is decomposed into a stack
103 : // of 4 Floating point arrays and the cleaning is done on all the the arrays
104 : // simultaneosly. The criterion for choosing the brightest pixel has been
105 : // generalised by using the "length" of the Stokesvector in 4 dimensional
106 : // space.
107 : //
108 : // A companion class to this one is MaskedClarkCleanLatModel. This provides
109 : // the same functionality but is used with MaskedArrays which indicate which
110 : // regions of the model to search for clean components.
111 : //
112 : // </synopsis>
113 : //
114 : // <example>
115 : // <srcblock>
116 : // casacore::Matrix<casacore::Float> psf(12,12), dirty(10,10), initialModel(10,10);
117 : // ...put appropriate values into psf, dirty, & initialModel....
118 : // ClarkCleanLatModel<casacore::Float> deconvolvedModel(initialModel);
119 : // ConvolutionEquation convEqn(psf, dirty);
120 : // deconvolvedModel.setGain(0.2);
121 : // deconvolvedModel.setNumberIterations(1000);
122 : // casacore::Bool convWorked = deconvolvedModel.solve(convEqn);
123 : // if (convWorked)
124 : // ConvEqn.residual(deconvolvedModel, finalResidual);
125 : // </srcblock>
126 : // </example>
127 : //
128 : // <motivation>
129 : // This class is needed to deconvolve images.
130 : // </motivation>
131 : //
132 : // <templating arg=T>
133 : // I have tested this class with Arrays of
134 : // <li> Float
135 : // <li> StokesVector
136 : // </templating>
137 : //
138 : // <todo asof="1996/05/02">
139 : // <li> Make changes so that multidimensions work as advertised
140 : // <li> compare timing with other clean implementations (ie, Mark's
141 : // CleanTools, SDE, AIPS & miriad)
142 : // </todo>
143 :
144 : class ClarkCleanLatModel:
145 : public LinearModel< casacore::Lattice<casacore::Float> >,
146 : public Iterate
147 : {
148 : public:
149 : // The default constructor does nothing more than initialise a zero length
150 : // array to hold the deconvolved model. If this constructor is used then
151 : // the actual model must be set using the setModel() function of the
152 : // LatticeModel class.
153 : ClarkCleanLatModel();
154 :
155 : // Construct the ClarkCleanLatModel object and initialise the model.
156 : ClarkCleanLatModel(casacore::Lattice<casacore::Float> & model);
157 :
158 : // Construct the ClarkCleanLatModel object and initialise the model ans mask
159 : ClarkCleanLatModel(casacore::Lattice<casacore::Float> & model, casacore::Lattice<casacore::Float> & mask);
160 :
161 : // Construct the ClarkCleanLatModel object and initialise the model ans mask
162 : ClarkCleanLatModel(casacore::Lattice<casacore::Float> & model, casacore::Lattice<casacore::Float> & residual,
163 : casacore::Lattice<casacore::Float> & mask);
164 : // Destroy!
165 : virtual ~ClarkCleanLatModel();
166 :
167 1982 : virtual const casacore::Lattice<casacore::Float> & getModel() const { return *itsModelPtr; }
168 : virtual void setModel(const casacore::Lattice<casacore::Float> & model);
169 : virtual void setModel(casacore::Lattice<casacore::Float> & model);
170 :
171 : const casacore::Lattice<casacore::Float> & getMask() const;
172 : void setMask(const casacore::Lattice<casacore::Float> & mask);
173 :
174 : void setResidual( casacore::Lattice<casacore::Float> & residual);
175 120 : virtual const casacore::Lattice<casacore::Float> & getResidual() const { return *itsResidualPtr; }
176 :
177 120 : casacore::Int getNumberIterations(){ return numberIterations(); }
178 :
179 120 : casacore::Float getMaxResidual() { return itsMaxRes;};
180 : // Using a Clark clean deconvolution proceedure solve for an improved
181 : // estimate of the deconvolved object. The convolution/residual equation
182 : // contains the psf and dirty image. When called with a ResidualEquation
183 : // arguement a quite general interface is used that is slow. The
184 : // convolution equation contains functions that speed things up. The
185 : // functions return false if the deconvolution could not be done.
186 : // <group>
187 : casacore::Bool solve(LatConvEquation & eqn);
188 : casacore::Bool singleSolve(LatConvEquation & eqn, casacore::Lattice<casacore::Float> & residual);
189 : // </group>
190 :
191 : // The user can be asked whether to stop after every minor cycle
192 : // <group>
193 : virtual void setChoose(const casacore::Bool askForChoice);
194 : virtual casacore::Bool getChoose();
195 : // </group>
196 :
197 : // These remaining functions set various "knobs" that the user can tweak and
198 : // are specific to the Clark clean algorithm. The more generic parameters
199 : // ie. clean gain, and maximum residual fluxlimit, are set using functions in
200 : // the Iterate base class. The get functions return the value that was
201 : // actually used after the cleaning was done.
202 :
203 : // set the size of the PSF used in the minor iterations. If not set it
204 : // defaults to the largest useful Psf (ie. min(modelSize*2, psfSize))
205 : // <group>
206 : virtual void setPsfPatchSize(const casacore::IPosition & psfPatchSize);
207 : virtual casacore::IPosition getPsfPatchSize();
208 : // </group>
209 :
210 : // Set the size of the histogram used to determine how many pixels are
211 : // "active" in a minor iteration. Default value of 1000 is OK for
212 : // everything except very small cleans.
213 : // <group>
214 : virtual void setHistLength(const casacore::uInt histBins);
215 : virtual casacore::uInt getHistLength();
216 : // </group>
217 :
218 : // Set the maximum number of minor iterations to perform for each major
219 : // cycle.
220 : // <group>
221 : virtual void setMaxNumberMinorIterations(const casacore::uInt maxNumMinorIterations);
222 : virtual casacore::uInt getMaxNumberMinorIterations();
223 : // </group>
224 :
225 : // Set and get the initial number of iterations
226 : // <group>
227 : virtual void setInitialNumberIterations(const casacore::uInt initialNumberIterations);
228 : virtual casacore::uInt getInitialNumberIterations();
229 : // </group>
230 :
231 : // Set the maximum number of major cycles to perform
232 : // <group>
233 : virtual void setMaxNumberMajorCycles(const casacore::uInt maxNumMajorCycles);
234 : virtual casacore::uInt getMaxNumberMajorCycles();
235 : // </group>
236 :
237 : // Set the maximum number of active pixels to use in the minor
238 : // iterations. The specified number can be exceeded if the topmost bin of
239 : // the histogram contains more pixels than specified here. The default is
240 : // 10,000 which is suitable for images of 512by512 pixels. Reduce this for
241 : // smaller images and increase it for larger ones.
242 : // <group>
243 : virtual void setMaxNumPix(const casacore::uInt maxNumPix );
244 : virtual casacore::uInt getMaxNumPix();
245 : // </group>
246 :
247 :
248 : // Set the maximum exterior psf value. This is used to determine when to
249 : // stop the minor itartions. This is normally determined from the psf and
250 : // the number set here is only used if this cannot be determined. The
251 : // default is zero.
252 : // <group>
253 : virtual void setMaxExtPsf(const casacore::Float maxExtPsf );
254 : virtual casacore::Float getMaxExtPsf();
255 : // </group>
256 :
257 : // The total flux density in the model.
258 : casacore::Float modelFlux();
259 :
260 : // An exponent on the F(m,n) factor (see Clark[1980]) which influences how
261 : // quickly active pixels are treated as unreliable. Larger values mean
262 : // more major iterations. The default is zero. I have no experience on
263 : // when to use this factor.
264 : // <group>
265 : virtual void setSpeedup(const casacore::Float speedup);
266 : virtual casacore::Float getSpeedup();
267 : // </group>
268 : //Set the cycle factor....the larger this is the shallower is the minor
269 : //cycle
270 : virtual void setCycleFactor(const casacore::Float factor);
271 :
272 :
273 : // Set/get the progress display
274 : // <group>
275 0 : virtual void setProgress(ClarkCleanProgress& ccp) { itsProgressPtr = &ccp; }
276 0 : virtual ClarkCleanProgress& getProgress() { return *itsProgressPtr; }
277 : // </group>
278 :
279 : private:
280 : // Do all the minor iterations for one major cycle. Cleaning stops
281 : // when the flux or iteration limit is reached.
282 : void doMinorIterations(CCList & activePixels,
283 : casacore::Matrix<casacore::Float> & psfPatch,
284 : casacore::Float fluxLimit,
285 : casacore::uInt & numberIterations,
286 : casacore::Float Fmn,
287 : const casacore::uInt totalIterations,
288 : casacore::Float& totalFlux);
289 :
290 : void cacheActivePixels(CCList & activePixels,
291 : const casacore::Lattice<casacore::Float> & residual, casacore::Float fluxLimit);
292 :
293 : // make histogram of absolute values in array
294 : void absHistogram(casacore::Vector<casacore::Int> & hist, casacore::Float & minVal,
295 : casacore::Float & maxVal, const casacore::Lattice<casacore::Float> & data);
296 :
297 : // Determine the flux limit if we only select the maxNumPix biggest
298 : // residuals. Flux limit is not exact due to quantising by the histogram
299 : casacore::Float biggestResiduals(casacore::Float & maxRes, const casacore::uInt maxNumPix,
300 : const casacore::Float fluxLimit,
301 : const casacore::Lattice<casacore::Float> & residual);
302 :
303 : // Work out the size of the Psf patch to use.
304 : casacore::Float getPsfPatch(casacore::Matrix<casacore::Float> & psfPatch, LatConvEquation & eqn);
305 :
306 : // The maximum residual is the absolute maximum.
307 : casacore::Float maxResidual(const casacore::Lattice<casacore::Float> & residual);
308 : void maxVect(casacore::Block<casacore::Float> & maxVal, casacore::Float & absVal, casacore::Int & offset,
309 : const CCList & activePixels);
310 : void subtractComponent(CCList & activePixels, const casacore::Block<casacore::Float> & maxVal,
311 : const casacore::Block<casacore::Int> & maxPos, const casacore::Matrix<casacore::Float> & psf);
312 : casacore::Float absMaxBeyondDist(const casacore::IPosition & maxDist, const casacore::IPosition & centre,
313 : const casacore::Lattice<casacore::Float> & psf);
314 : casacore::Bool stopnow();
315 : casacore::Int getbig(casacore::Float const * pixValPtr, casacore::Int const * pixPosPtr, const casacore::Int nPix,
316 : const casacore::Float fluxLimit,
317 : const casacore::Float * const residualPtr, const casacore::Float * const maskPtr,
318 : const casacore::uInt npol, const casacore::Int nx, const casacore::Int ny);
319 :
320 : void updateModel(CCList & cleanComponents);
321 :
322 : casacore::Lattice<casacore::Float> * itsModelPtr;
323 : casacore::Lattice<casacore::Float> * itsResidualPtr;
324 : const casacore::Lattice<casacore::Float> * itsSoftMaskPtr;
325 : casacore::uInt itsMaxNumPix;
326 : casacore::uInt itsHistBins;
327 : casacore::uInt itsMaxNumberMinorIterations;
328 : casacore::uInt itsInitialNumberIterations;
329 : casacore::Int itsMaxNumberMajorCycles;
330 : casacore::Float itsMaxExtPsf;
331 : casacore::Float itsMaxRes;
332 : casacore::IPosition itsPsfPatchSize;
333 : casacore::Bool itsChoose;
334 : casacore::Float itsSpeedup;
335 : casacore::Float itsCycleFactor;
336 : casacore::LogIO itsLog;
337 : ClarkCleanProgress* itsProgressPtr;
338 : casacore::Bool itsJustStarting;
339 : casacore::Bool itsWarnFlag;
340 :
341 : casacore::Bool itsLocalResTL;
342 :
343 : };
344 :
345 :
346 : } //# NAMESPACE CASA - END
347 :
348 : #endif
|
