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