Line data Source code
1 : //# Cleaner.h: this defines Cleaner a class for doing convolution
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,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: casa-feedback@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: LatticeCleaner.h 20739 2009-09-29 01:15:15Z Malte.Marquarding $
28 :
29 : #ifndef SYNTHESIS_MATRIXCLEANER_H
30 : #define SYNTHESIS_MATRIXCLEANER_H
31 :
32 : //# Includes
33 : #include <casacore/casa/aips.h>
34 : #include <casacore/casa/Quanta/Quantum.h>
35 : #include <casacore/lattices/Lattices/TempLattice.h>
36 : #include <casacore/casa/Arrays/IPosition.h>
37 : #include <casacore/casa/Arrays/Vector.h>
38 : #include <casacore/casa/Containers/Block.h>
39 : #include <casacore/lattices/LatticeMath/LatticeCleaner.h>
40 : #include <casacore/casa/Arrays/ArrayFwd.h>
41 : #include <casacore/casa/Utilities/CountedPtr.h>
42 :
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 :
45 : //# Forward Declarations
46 : class MatrixNACleaner;
47 : // <summary>A copy of casacore::LatticeCleaner but just using 2-D matrices</summary>
48 : // <synopsis> It is a mimic of the casacore::LatticeCleaner class but avoid a lot of
49 : // of the lattice to array and back copies and uses openmp in the obvious places
50 : // </synopsis>
51 :
52 : // <summary>A class for doing multi-dimensional cleaning</summary>
53 :
54 : // <use visibility=export>
55 :
56 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="tLatticeCleaner">
57 : // </reviewed>
58 :
59 : // <prerequisite>
60 : // <li> The mathematical concept of deconvolution
61 : // </prerequisite>
62 : //
63 : // <etymology>
64 :
65 : // The MatrixCleaner class will deconvolve 2-D arrays of floats.
66 :
67 : // </etymology>
68 : //
69 : // <synopsis>
70 : // This class will perform various types of Clean deconvolution
71 : // on Lattices.
72 : //
73 : // </synopsis>
74 : //
75 : // <example>
76 : // <srcblock>
77 : // </srcblock>
78 : // </example>
79 : //
80 : // <motivation>
81 : // </motivation>
82 : //
83 : // <thrown>
84 : // <li> casacore::AipsError: if psf has more dimensions than the model.
85 : // </thrown>
86 : //
87 : // <todo asof="yyyy/mm/dd">
88 : // </todo>
89 :
90 : class MatrixCleaner
91 : {
92 : public:
93 :
94 : // Create a cleaner : default constructor
95 : MatrixCleaner();
96 :
97 : // Create a cleaner for a specific dirty image and PSF
98 : MatrixCleaner(const casacore::Matrix<casacore::Float> & psf, const casacore::Matrix<casacore::Float> & dirty);
99 :
100 : // The copy constructor uses reference semantics
101 : MatrixCleaner(const MatrixCleaner& other);
102 :
103 : // The assignment operator also uses reference semantics
104 : MatrixCleaner & operator=(const MatrixCleaner & other);
105 :
106 : // The destructor does nothing special.
107 : ~MatrixCleaner();
108 :
109 : //just define the scales...nothing else is done
110 : //the user will need to call setPsf+makePsfScales+setDirty+makeDirtyScales
111 : //to be in a good state to clean.
112 : void defineScales(const casacore::Vector<casacore::Float>& scales);
113 :
114 : //Set the dirty image without calculating convolutions..
115 : //can be done by calling makeDirtyScales or setscales if one want to redo the
116 : //psfscales too.
117 : void setDirty(const casacore::Matrix<casacore::Float>& dirty);
118 : //Calculate the convolutions for the dirt
119 : //Obviously the
120 : void makeDirtyScales();
121 : // Update the dirty image only (equiv of setDirty + makeDirtyScales)
122 : void update(const casacore::Matrix<casacore::Float> & dirty);
123 :
124 : //change the psf
125 : //don't forget to redo the setscales or run makePsfScales,
126 : //followed by makeDirtyScales
127 : void setPsf(const casacore::Matrix<casacore::Float>& psf);
128 : //calculate the convolutions of the psf
129 : void makePsfScales();
130 :
131 : // Set a number of scale sizes. The units of the scale are pixels.
132 : // The 2 functions below assume you have the dirty image and the psf set
133 : // the convolutions are calculated automatically and the masks ones too
134 : // if it is set ....kept so as to be compatible function wise with LatticeCleaner
135 : casacore::Bool setscales(const casacore::Int nscales, const casacore::Float scaleInc=1.0);
136 :
137 : // Set a specific set of scales
138 : casacore::Bool setscales(const casacore::Vector<casacore::Float> & scales);
139 :
140 :
141 :
142 : // Set up control parameters
143 : // cleanType - type of the cleaning algorithm to use (HOGBOM, MULTISCALE)
144 : // niter - number of iterations
145 : // gain - loop gain used in cleaning (a fraction of the maximum
146 : // subtracted at every iteration)
147 : // aThreshold - absolute threshold to stop iterations
148 : // fThreshold - fractional threshold (i.e. given w.r.t. maximum residual)
149 : // to stop iterations. This parameter is specified as
150 : // casacore::Quantity so it can be given in per cents.
151 : // choose - unused at the moment, specify false. Original meaning is
152 : // to allow interactive decision on whether to continue iterations.
153 : // This method always returns true.
154 : casacore::Bool setcontrol(casacore::CleanEnums::CleanType cleanType, const casacore::Int niter,
155 : const casacore::Float gain, const casacore::Quantity& aThreshold,
156 : const casacore::Quantity& fThreshold);
157 :
158 : // This version of the method disables stopping on fractional threshold
159 : casacore::Bool setcontrol(casacore::CleanEnums::CleanType cleanType, const casacore::Int niter,
160 : const casacore::Float gain, const casacore::Quantity& threshold);
161 :
162 : // return how many iterations we did do
163 0 : casacore::Int iteration() const { return itsIteration; }
164 0 : casacore::Int numberIterations() const { return itsIteration; }
165 :
166 : // what iteration number to start on
167 0 : void startingIteration(const casacore::Int starting = 0) {itsStartingIter = starting; }
168 :
169 : //Total flux accumulated so far
170 : casacore::Float totalFlux() const {return itsTotalFlux;}
171 :
172 :
173 : // Clean an image.
174 : //return value gives you a hint of what's happening
175 : // 1 = converged
176 : // 0 = not converged but behaving normally
177 : // -1 = not converged and stopped on cleaning consecutive smallest scale
178 : // -2 = not converged and either large scale hit negative or diverging
179 : // -3 = clean is diverging rather than converging
180 : casacore::Int clean(casacore::Matrix<casacore::Float> & model, casacore::Bool doPlotProgress=false);
181 :
182 : // Set the mask
183 : // mask - input mask lattice
184 : // maskThreshold - if positive, the value is treated as a threshold value to determine
185 : // whether a pixel is good (mask value is greater than the threshold) or has to be
186 : // masked (mask value is below the threshold). Negative threshold switches mask clipping
187 : // off. The mask value is used to weight the flux during cleaning. This mode is used
188 : // to implement cleaning based on the signal-to-noise as opposed to the standard cleaning
189 : // based on the flux. The default threshold value is 0.9, which ensures the behavior of the
190 : // code is exactly the same as before this parameter has been introduced.
191 0 : void setMask(casacore::Matrix<casacore::Float> & mask, const casacore::Float& maskThreshold = 0.9);
192 : // Call the function below if the psf is changed ..no need to setMask again
193 : casacore::Bool makeScaleMasks();
194 :
195 : // remove the mask;
196 : // useful when keeping object and sending a new dirty image to clean
197 : // one can set another mask then
198 : void unsetMask();
199 :
200 : // Tell the algorithm to NOT clean just the inner quarter
201 : // (This is useful when multiscale clean is being used
202 : // inside a major cycle for MF or WF algorithms)
203 : // if true, the full image deconvolution will be attempted
204 0 : void ignoreCenterBox(casacore::Bool huh) { itsIgnoreCenterBox = huh; }
205 :
206 : // Consider the case of a point source:
207 : // the flux on all scales is the same, and the first scale will be chosen.
208 : // Now, consider the case of a point source with a *little* bit of extended structure:
209 : // thats right, the largest scale will be chosen. In this case, we should provide some
210 : // bias towards the small scales, or against the large scales. We do this in
211 : // an ad hoc manner, multiplying the maxima found at each scale by
212 : // 1.0 - itsSmallScaleBias * itsScaleSizes(scale)/itsScaleSizes(nScalesToClean-1);
213 : // Typical bias values range from 0.2 to 1.0.
214 0 : void setSmallScaleBias(const casacore::Float x=0.5) { itsSmallScaleBias = x; }
215 :
216 : // During early iterations of a cycled casacore::MS Clean in mosaicing, it common
217 : // to come across an ocsilatory pattern going between positive and
218 : // negative in the large scale. If this is set, we stop at the first
219 : // negative in the largest scale.
220 0 : void stopAtLargeScaleNegative() {itsStopAtLargeScaleNegative = true; }
221 :
222 : // Some algorithms require that the cycles be terminated when the image
223 : // is dominated by point sources; if we get nStopPointMode of the
224 : // smallest scale components in a row, we terminate the cycles
225 0 : void stopPointMode(casacore::Int nStopPointMode) {itsStopPointMode = nStopPointMode; }
226 :
227 : // After completion of cycle, querry this to find out if we stopped because
228 : // of stopPointMode
229 0 : casacore::Bool queryStopPointMode() const {return itsDidStopPointMode; }
230 :
231 : // speedup() will speed the clean iteration by raising the
232 : // threshold. This may be required if the threshold is
233 : // accidentally set too low (ie, lower than can be achieved
234 : // given errors in the approximate PSF).
235 : //
236 : // threshold(iteration) = threshold(0)
237 : // * ( exp( (iteration - startingiteration)/Ndouble )/ 2.718 )
238 : // If speedup() is NOT invoked, no effect on threshold
239 : void speedup(const casacore::Float Ndouble);
240 :
241 : // Look at what WE think the residuals look like
242 : // Assumes the first scale is zero-sized
243 : casacore::Matrix<casacore::Float> residual() { return itsDirtyConvScales[0]; }
244 : //slightly better approximation of the residual: it convolves the given model
245 : //with the psf and remove it from the dirty image put in setdirty
246 : casacore::Matrix<casacore::Float> residual(const casacore::Matrix<casacore::Float>& model);
247 : // Method to return threshold, including any speedup factors
248 : casacore::Float threshold() const;
249 :
250 : // Method to return the strength optimum achieved at the last clean iteration
251 : // The output of this method makes sense only if it is called after clean
252 0 : casacore::Float strengthOptimum() const { return itsStrengthOptimum; }
253 :
254 : // Helper function to optimize adding
255 : //static void addTo(casacore::Matrix<casacore::Float>& to, const casacore::Matrix<casacore::Float> & add);
256 :
257 : protected:
258 : friend class MatrixNACleaner;
259 : // Make sure that the peak of the Psf is within the image
260 : casacore::Bool validatePsf(const casacore::Matrix<casacore::Float> & psf);
261 :
262 : // Make an array of the specified scale
263 : void makeScale(casacore::Matrix<casacore::Float>& scale, const casacore::Float& scaleSize);
264 :
265 : // Make Spheroidal function for scale images
266 : casacore::Float spheroidal(casacore::Float nu);
267 :
268 : // Find the Peak of the matrix
269 : static casacore::Bool findMaxAbs(const casacore::Matrix<casacore::Float>& lattice,
270 : casacore::Float& maxAbs, casacore::IPosition& posMax);
271 :
272 : // This is made static since findMaxAbs is static(!).
273 : // Why is findMaxAbs static???
274 : // --SB
275 : static casacore::Bool findPSFMaxAbs(const casacore::Matrix<casacore::Float>& lattice,
276 : casacore::Float& maxAbs, casacore::IPosition& posMax,
277 : const casacore::Int& supportSize=100);
278 :
279 : casacore::Int findBeamPatch(const casacore::Float maxScaleSize, const casacore::Int& nx, const casacore::Int& ny,
280 : const casacore::Float psfBeam=4.0, const casacore::Float nBeams=20.0);
281 : // Find the Peak of the lattice, applying a mask
282 : casacore::Bool findMaxAbsMask(const casacore::Matrix<casacore::Float>& lattice, const casacore::Matrix<casacore::Float>& mask,
283 : casacore::Float& maxAbs, casacore::IPosition& posMax);
284 :
285 : // Helper function to reduce the box sizes until the have the same
286 : // size keeping the centers intact
287 : static void makeBoxesSameSize(casacore::IPosition& blc1, casacore::IPosition& trc1,
288 : casacore::IPosition &blc2, casacore::IPosition& trc2);
289 :
290 :
291 : casacore::CleanEnums::CleanType itsCleanType;
292 : casacore::Float itsGain;
293 : casacore::Int itsMaxNiter; // maximum possible number of iterations
294 : casacore::Quantum<casacore::Double> itsThreshold;
295 : casacore::CountedPtr<casacore::Matrix<casacore::Float> > itsMask;
296 : casacore::IPosition itsPositionPeakPsf;
297 : casacore::Float itsSmallScaleBias;
298 : casacore::Block<casacore::Matrix<casacore::Float> > itsScaleMasks;
299 : casacore::Block<casacore::Matrix<casacore::Complex> > itsScaleXfrs;
300 : casacore::Bool itsScalesValid;
301 : casacore::Int itsNscales;
302 : casacore::Float itsMaskThreshold;
303 :
304 : //# The following functions are used in various places in the code and are
305 : //# documented in the .cc file. Static functions are used when the functions
306 : //# do not modify the object state. They ensure that implicit assumptions
307 : //# about the current state and implicit side-effects are not possible
308 : //# because all information must be supplied in the input arguments
309 :
310 :
311 : casacore::CountedPtr<casacore::Matrix<casacore::Float> > itsDirty;
312 : casacore::CountedPtr<casacore::Matrix<casacore::Complex> >itsXfr;
313 :
314 : casacore::Vector<casacore::Float> itsScaleSizes;
315 :
316 : casacore::Block<casacore::Matrix<casacore::Float> > itsScales;
317 : casacore::Block<casacore::Matrix<casacore::Float> > itsPsfConvScales;
318 : casacore::Block<casacore::Matrix<casacore::Float> > itsDirtyConvScales;
319 :
320 :
321 : casacore::Int itsIteration; // what iteration did we get to?
322 : casacore::Int itsStartingIter; // what iteration did we get to?
323 : casacore::Quantum<casacore::Double> itsFracThreshold;
324 :
325 : casacore::Float itsMaximumResidual;
326 : casacore::Float itsStrengthOptimum;
327 :
328 :
329 : casacore::Vector<casacore::Float> itsTotalFluxScale;
330 : casacore::Float itsTotalFlux;
331 :
332 : // Calculate index into PsfConvScales
333 : casacore::Int index(const casacore::Int scale, const casacore::Int otherscale);
334 :
335 : casacore::Bool destroyScales();
336 : casacore::Bool destroyMasks();
337 :
338 : casacore::Bool itsIgnoreCenterBox;
339 : casacore::Bool itsStopAtLargeScaleNegative;
340 : casacore::Int itsStopPointMode;
341 : casacore::Bool itsDidStopPointMode;
342 :
343 : casacore::IPosition psfShape_p;
344 : casacore::Bool noClean_p;
345 :
346 : private:
347 :
348 : // casacore::Memory to be allocated per TempLattice
349 : casacore::Double itsMemoryMB;
350 :
351 : // Let the user choose whether to stop
352 : casacore::Bool itsChoose;
353 :
354 : // Threshold speedup factors:
355 : casacore::Bool itsDoSpeedup; // if false, threshold does not change with iteration
356 : casacore::Float itsNDouble;
357 :
358 : //# Stop now?
359 : //#// casacore::Bool stopnow(); Removed on 8-Apr-2004 by GvD
360 :
361 : casacore::Bool itsJustStarting;
362 :
363 : // threshold for masks. If negative, mask values are used as weights and no pixels are
364 : // discarded (although effectively they would be discarded if the mask value is 0.)
365 : };
366 :
367 : } //# NAMESPACE CASA - END
368 :
369 : #endif
|
