Line data Source code
1 : //# ClarkCleanModel.h: this defines ClarkCleanModel
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_CLARKCLEANMODEL_H
30 : #define SYNTHESIS_CLARKCLEANMODEL_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <synthesis/MeasurementEquations/ArrayModel.h>
37 : #include <synthesis/MeasurementEquations/Iterate.h>
38 : //#include <synthesis/MeasurementEquations/ResidualEquation.h>
39 : #include <synthesis/MeasurementEquations/ConvolutionEquation.h>
40 : #include <casacore/casa/Logging/LogIO.h>
41 :
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 :
44 : class ClarkCleanProgress;
45 :
46 : // <summary>
47 : // A Class for performing the Clark Clean Algorithm on Arrays
48 : // </summary>
49 :
50 : // <use visibility=export>
51 :
52 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
53 : // </reviewed>
54 :
55 : // <prerequisite>
56 : // <li> ResidualEquation/ConvolutionEquation
57 : // <li> LinearModel/LinearEquation Paradigm
58 : // </prerequisite>
59 : //
60 : // <etymology>
61 : // This class is called ClarkCleanModel because thats the algorithm it uses
62 : // deconvolve the model.
63 : // </etymology>
64 : //
65 : // <synopsis>
66 : // This class is used to perform the Clark Clean Algorithm on an
67 : // Array. Only the deconvolved model of the sky are directly stored by this
68 : // class. The point spread function (psf) and convolved (dirty) image are
69 : // stored in a companion class which is must be derived from
70 : // ResidualEquation.
71 : //
72 : // The cleaning works like this. The user constructs a ClarkCleanModel by
73 : // specifying an initial model of the sky. This can by be
74 : // one,two,three... dimensional depending on the dimension of the psf (see
75 : // below). The user then constructs a class which implements the forward
76 : // equation between the model and the dirty image. Typically this will be
77 : // the ConvolutionEquation class, although any class which has a
78 : // ResidualEquation interface will be work (but perhaps very slowly, as the
79 : // ConvolutionEquation class has member functions optimised for cleaning)
80 : //
81 : // The user then calls the solve() function (with the appropriate equation
82 : // class as an arguement), and this class will perform the Clark clean.
83 : // The various clean parameters are set (prior to calling solve) using the
84 : // functions derived from the Iterate class, in particular setGain(),
85 : // setNumberIterations() & setThreshold() (to set a flux limit).
86 : //
87 : // The solve() function does not return either the deconvolved model or the
88 : // residuals. The solved model can be obtained using the getModel() function
89 : // (derived from ArrayModel()) and the residual can be obtained using the
90 : // residual() member function of the Convolution/Residual Equation Class.
91 : //
92 : // The size and shape of the model used in this class MUST be the same as
93 : // the convolved data (Dirty Image), stored in the companion
94 : // ResidualEquation Class. However the model (and convolved data) can have
95 : // more dimensions than the psf, as well as a different size (either larger
96 : // or smaller). When the dimensionality is different the cleaning is done
97 : // independendtly in each "plane" of the model. (Note this has not
98 : // been implemented yet but is relatively simple to do if necessary).
99 : //
100 : // This multi-dimensionalty is exploited when cleaning arrays of
101 : // StokesVectors. Here the casacore::Array of StokesVectors is decomposed into a stack
102 : // of 4 Floating point arrays and the cleaning is done on all the the arrays
103 : // simultaneosly. The criterion for choosing the brightest pixel has been
104 : // generalised by using the "length" of the Stokesvector in 4 dimensional
105 : // space.
106 : //
107 : // A companion class to this one is MaskedClarkCleanModel. This provides
108 : // the same functionality but is used with MaskedArrays which indicate which
109 : // regions of the model to search for clean components.
110 : //
111 : // </synopsis>
112 : //
113 : // <example>
114 : // <srcblock>
115 : // casacore::Matrix<casacore::Float> psf(12,12), dirty(10,10), initialModel(10,10);
116 : // ...put appropriate values into psf, dirty, & initialModel....
117 : // ClarkCleanModel<casacore::Float> deconvolvedModel(initialModel);
118 : // ConvolutionEquation convEqn(psf, dirty);
119 : // deconvolvedModel.setGain(0.2);
120 : // deconvolvedModel.setNumberIterations(1000);
121 : // casacore::Bool convWorked = deconvolvedModel.solve(convEqn);
122 : // casacore::Array<casacore::Float> finalModel, residuals;
123 : // if (convWorked){
124 : // finalModel = deconvolvedModel.getModel();
125 : // ConvEqn.residual(deconvolvedModel, finalResidual);
126 : // }
127 : // </srcblock>
128 : // </example>
129 : //
130 : // <motivation>
131 : // This class is needed to deconvolve images.
132 : // </motivation>
133 : //
134 : // <templating arg=T>
135 : // I have tested this class with Arrays of
136 : // <li> Float
137 : // <li> StokesVector
138 : // </templating>
139 : //
140 : // <todo asof="1996/05/02">
141 : // <li> Make changes so that multidimensions work as advertised
142 : // <li> compare timing with other clean implementations (ie, Mark's
143 : // CleanTools, SDE, AIPS & miriad)
144 : // </todo>
145 :
146 : class ClarkCleanModel:
147 : public ArrayModel<casacore::Float>,
148 : public Iterate
149 : {
150 : public:
151 : // The default constructor does nothing more than initialise a zero length
152 : // array to hold the deconvolved model. If this constructor is used then
153 : // the actual model must be set using the setModel() function of the
154 : // ArrayModel class.
155 : ClarkCleanModel();
156 : // Construct the ClarkCleanModel object and initialise the model.
157 : ClarkCleanModel(casacore::Array<casacore::Float> & model);
158 : // Construct the ClarkCleanModel object and initialise the model ans mask
159 : ClarkCleanModel(casacore::Array<casacore::Float> & model, casacore::Array<casacore::Float> & mask);
160 :
161 : void getMask(casacore::Array<casacore::Float>& mask) const;
162 : void setMask(const casacore::Array<casacore::Float>& mask);
163 : void setMask(casacore::Array<casacore::Float> & mask);
164 :
165 : void getModel(casacore::Array<casacore::Float>& model) const;
166 : void setModel(const casacore::Array<casacore::Float>& model);
167 : void setModel(casacore::Array<casacore::Float> & model);
168 :
169 : // Set/get the progress display
170 : // <group>
171 0 : virtual void setProgress(ClarkCleanProgress& ccp) { itsProgressPtr = &ccp; }
172 0 : virtual ClarkCleanProgress& getProgress() { return *itsProgressPtr; }
173 : // </group>
174 :
175 : // Using a Clark clean deconvolution proceedure solve for an improved
176 : // estimate of the deconvolved object. The convolution/residual equation
177 : // contains the psf and dirty image. When called with a ResidualEquation
178 : // arguement a quite general interface is used that is slow. The
179 : // convolution equation contains functions that speed things up. The
180 : // functions return false if the deconvolution could not be done.
181 : // <group>
182 : casacore::Bool solve(ConvolutionEquation & eqn);
183 : casacore::Bool singleSolve(ConvolutionEquation & eqn, casacore::Array<casacore::Float>& residual);
184 : // </group>
185 :
186 : // These functions set various "knobs" that the user can tweak and are
187 : // specific to the Clark clean algorithm. The more generic parameters
188 : // ie. clean gain, and maximum residual fluxlimit, are set using functions
189 : // in the Iterate base class. The get functions return the value that was
190 : // actually used after the cleaning was done.
191 : // <group>
192 : // set the size of the PSF used in the minor iterations. If not set it
193 : // defaults to the largest useful Psf (ie. min(modelSize*2, psfSize))
194 : virtual void setPsfPatchSize(const casacore::IPosition & psfPatchSize);
195 : virtual casacore::IPosition getPsfPatchSize();
196 : // Set the size of the histogram used to determine how many pixels are
197 : // "active" in a minor iteration. Default value is 1000 is OK for
198 : // everything except very small cleans.
199 : virtual void setHistLength(const casacore::uInt HistBins );
200 : virtual casacore::uInt getHistLength();
201 : // Set the maximum number of minor iterations to perform for each major
202 : // cycle.
203 : virtual void setMaxNumberMinorIterations(const casacore::uInt maxNumMinorIterations);
204 : virtual casacore::uInt getMaxNumberMinorIterations();
205 : // Set and get the initial number of iterations
206 : virtual void setInitialNumberIterations(const casacore::uInt initialNumberIterations);
207 : virtual casacore::uInt getInitialNumberIterations();
208 : // Set the maximum number of major cycles to perform
209 : virtual void setMaxNumberMajorCycles(const casacore::uInt maxNumMajorCycles);
210 : virtual casacore::uInt getMaxNumberMajorCycles();
211 : // Set the maximum number of active pixels to use in the minor
212 : // iterations. The specified number can be exceeded if the topmost bin of
213 : // the histogram contains more pixels than specified here. The default is
214 : // 10,000 which is suitable for images of 512by512 pixels. Reduce this for
215 : // smaller images and increase it for larger ones.
216 : virtual void setMaxNumPix(const casacore::uInt maxNumPix );
217 : virtual casacore::uInt getMaxNumPix();
218 : // Set the maximum exterior psf value. This is used to determine when to
219 : // stop the minor itartions. This is normally determined from the psf and
220 : // the number set here is only used if this cannot be determined. The
221 : // default is zero.
222 : virtual void setMaxExtPsf(const casacore::Float maxExtPsf );
223 : virtual casacore::Float getMaxExtPsf();
224 : // An exponent on the F(m,n) factor (see Clark[1980]) which influences how
225 : // quickly active pixels are treated as unreliable. Larger values mean
226 : // more major iterations. The default is zero. I have no experience on
227 : // when to use this factor.
228 : virtual void setSpeedup(const casacore::Float speedup );
229 : virtual casacore::Float getSpeedup();
230 : // The structure of various AIPS++ algorithms creates major cycles around
231 : // the Clark Clean (or other deconvolution algrithms. The cycleSpeedup
232 : // parameter causes the threshold to edge up as
233 : // thresh = thresh_0 * 2^( iter/cycleSpeedup );
234 : // ignored if cycleSpeedup <= 0.
235 : virtual void setCycleSpeedup(const casacore::Float speedup );
236 : virtual casacore::Float getCycleSpeedup();
237 : // We are overwriting Iterate's threshold() method to put out speedup in it
238 : virtual casacore::Float threshold();
239 : // The user can be asked whether to stop after every minor cycle
240 : virtual void setChoose(const casacore::Bool askForChoice);
241 : virtual casacore::Bool getChoose();
242 : // </group>
243 0 : casacore::Float getMaxResidual() { return itsMaxRes;};
244 :
245 : private:
246 : // Do all the minor iterations for one major cycle. Cleaning stops
247 : // when the flux or iteration limit is reached.
248 : void doMinorIterations(casacore::Array<casacore::Float> & model,
249 : casacore::Matrix<casacore::Float> & pixelValue,
250 : const casacore::Matrix<casacore::Int> & pixelPos,
251 : const casacore::Int numPix,
252 : casacore::Matrix<casacore::Float> & psfPatch,
253 : casacore::Float fluxLimit,
254 : casacore::uInt & numberIterations,
255 : casacore::Float Fmn,
256 : const casacore::uInt totalIterations,
257 : casacore::Float& totalflux);
258 : // Find all the pixels in the residual that are greater than fluxlimit, and
259 : // store the values in the pixelsValue casacore::Matrix, and their positions in the
260 : // pixelPos Matrix. Increases the size of the output matrices as
261 : // necessary, but does not decrease them. So the actual number of "active"
262 : // pixels is returned. This will always be less than (or equal to) the matrix
263 : // size.
264 : casacore::Int cacheActivePixels(casacore::Matrix<casacore::Float> & pixVal, casacore::Matrix<casacore::Int> & pixPos,
265 : const casacore::Array<casacore::Float> & data, const casacore::Float fluxLimit);
266 : // make histogram of absolute values in array
267 : void absHistogram(casacore::Vector<casacore::Int> & hist, casacore::Float & minVal,
268 : casacore::Float & maxVal, const casacore::Array<casacore::Float> & data);
269 : // Determine the flux limit if we only select the maxNumPix biggest
270 : // residuals. Flux limit is not exact due to quantising by the histogram
271 : casacore::Float biggestResiduals(casacore::Float & maxRes, const casacore::uInt maxNumPix,
272 : const casacore::Float fluxLimit, const casacore::Array<casacore::Float> & residual);
273 : // Work out the size of the Psf patch to use.
274 : casacore::Float getPsfPatch(casacore::Array<casacore::Float>& psfPatch, ConvolutionEquation& eqn);
275 : // The maximum residual is the absolute maximum.
276 : casacore::Float maxResidual(const casacore::Array<casacore::Float> & residual);
277 : void maxVect(casacore::Vector<casacore::Float> & maxVal, casacore::Float & absVal, casacore::Int & offset,
278 : const casacore::Matrix<casacore::Float> & pixVal, const casacore::Int numPix);
279 : void subtractComponent(casacore::Matrix<casacore::Float> & pixVal, const casacore::Matrix<casacore::Int> & pixPos,
280 : const casacore::Int numPix, const casacore::Vector<casacore::Float> & maxVal,
281 : const casacore::Vector<casacore::Int> & maxPos, const casacore::Matrix<casacore::Float> & psf);
282 : casacore::Float absMaxBeyondDist(const casacore::IPosition &maxDist, const casacore::IPosition ¢re,
283 : const casacore::Array<casacore::Float> &array);
284 : casacore::Bool stopnow();
285 :
286 : casacore::uInt theHistBins;
287 : casacore::Float theMaxExtPsf;
288 : casacore::uInt theMaxNumberMinorIterations;
289 : casacore::uInt theInitialNumberIterations;
290 : casacore::Int theMaxNumberMajorCycles;
291 : casacore::uInt theMaxNumPix;
292 : casacore::IPosition thePsfPatchSize;
293 : casacore::Float theSpeedup;
294 : casacore::Float theCycleSpeedup;
295 : casacore::Bool theChoose;
296 : casacore::Array<casacore::Float> theMask;
297 : casacore::LogIO theLog;
298 : // There are too many iterations counters.
299 : // This one is required for theCycleSpeedup and threshold(),
300 : // and just keeps track of the total iterations done by THIS
301 : // ClarkCleanModel
302 : casacore::Int theIterCounter;
303 : ClarkCleanProgress* itsProgressPtr;
304 : casacore::Bool itsJustStarting;
305 :
306 : casacore::Float itsMaxRes;
307 : };
308 :
309 :
310 : } //# NAMESPACE CASA - END
311 :
312 : #endif
|
