Line data Source code
1 : //# IncCEMemModel.h: this defines IncCEMemModel
2 : //# Copyright (C) 1996,1997,1998,1999,2000
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 : //# $Id$
27 :
28 : #ifndef SYNTHESIS_INCCEMEMMODEL_H
29 : #define SYNTHESIS_INCCEMEMMODEL_H
30 :
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 <casacore/lattices/Lattices/Lattice.h>
37 : #include <casacore/images/Images/PagedImage.h>
38 : #include <synthesis/MeasurementEquations/IncEntropy.h>
39 : #include <synthesis/MeasurementEquations/LinearEquation.h>
40 : #include <synthesis/MeasurementEquations/LinearModel.h>
41 : #include <casacore/casa/Logging/LogIO.h>
42 :
43 :
44 : namespace casa { //# NAMESPACE CASA - BEGIN
45 :
46 : // Forward declaration
47 : class LatConvEquation;
48 : class CEMemProgress;
49 :
50 :
51 : // <summary> performs MEM algorithm incrementally </summary>
52 :
53 : // <use visibility=export>
54 :
55 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
56 : // </reviewed>
57 :
58 : // <prerequisite>
59 : // <li> ResidualEquation/ConvolutionEquation
60 : // <li> LinearModel/LinearEquation Paradigm
61 : // </prerequisite>
62 : //
63 : // <etymology>
64 : // This class is called CEMemModel because it uses the Cornwell and
65 : // Evans MEM algorithm to deconvolve the model. The "Inc" is from Incremental,
66 : // as the gradient is calculated from an incremental dirty image, but the
67 : // entropy is calculated from the previous model plus the deltaModel
68 : // </etymology>
69 : //
70 : // <synopsis>
71 : // This class is used to perform the Cornwell and Evans MEM Algorithm on an
72 : // Array. Only the deconvolved model of the sky are directly stored by this
73 : // class. The point spread function (psf) and convolved (dirty) image are
74 : // stored in a companion class which is must be derived from
75 : // ResidualEquation.
76 : //
77 : // The deconvolution works like this. The user constructs a CEMemModel by
78 : // specifying an initial model of the sky. This can by be
79 : // one,two,three... dimensional depending on the dimension of the psf (see
80 : // below). The user then constructs a class which implements the forward
81 : // equation between the model and the dirty image. Typically this will be
82 : // the ConvolutionEquation class, although any class which has a
83 : // ResidualEquation interface will work (but perhaps very slowly, as the
84 : // ConvolutionEquation class has member functions optimised for CLEAN and MEM)
85 : //
86 : // The user then calls the solve() function (with the appropriate equation
87 : // class as an arguement), and this class will perform the MEM deconvolution.
88 : // The various MEM parameters are set (prior to calling solve) using the
89 : // functions derived from the Iterate class, in particular
90 : // setNumberIterations().
91 : //
92 : // The solve() function does not return either the deconvolved model or the
93 : // residuals. The solved model can be obtained using the getModel() function
94 : // (derived from ArrayModel()) and the residual can be obtained using the
95 : // residual() member function of the Convolution/Residual Equation Class.
96 : //
97 : // The size and shape of the model used in this class MUST be the same as
98 : // the convolved data (Dirty Image), stored in the companion
99 : // ResidualEquation Class. However the model (and convolved data) can have
100 : // more dimensions than the psf, as well as a different size (either larger
101 : // or smaller). When the dimensionality is different the deconvolution is done
102 : // independendtly in each "plane" of the model. (Note this has not
103 : // been implemented yet but is relatively simple to do if necessary).
104 : //
105 : // StokesVectors are not yet implemented.
106 : //
107 : // A companion class to this one is MaskedCEMemModel. This provides
108 : // the same functionality but is used with MaskedArrays which indicate which
109 : // regions of the NewtonRaphson (residual) image to apply when forming the
110 : // step image (MaskedCEMemModel is not yet implemented).
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 : // CEMemModel<casacore::Float> deconvolvedModel(initialModel);
119 : // ConvolutionEquation convEqn(psf, dirty);
120 : // deconvolvedModel.setSigma(0.001);
121 : // deconvolvedModel.setTargetFlux(-2.500);
122 : // deconvolvedModel.setNumberIterations(30);
123 : // casacore::Bool convWorked = deconvolvedModel.solve(convEqn);
124 : // casacore::Array<casacore::Float> finalModel, residuals;
125 : // if (convWorked){
126 : // finalModel = deconvolvedModel.getModel();
127 : // ConvEqn.residual(deconvolvedModel, finalResidual);
128 : // }
129 : // </srcblock>
130 : // </example>
131 : //
132 : // <motivation>
133 : // This class is needed to deconvolve extended images.
134 : // In time, the MEM algorithm will be a principle player in the
135 : // mosaicing stuff.
136 : // </motivation>
137 : //
138 : // <templating arg=T>
139 : // For testing:
140 : // <li> casacore::Float: lets try it first
141 : // <li> StokesVector: will require lots more work
142 : // </templating>
143 : //
144 : // <todo asof="1998/12/02">
145 : // <li> We need to implement soft Masking
146 : // </todo>
147 :
148 : class IncCEMemModel: public LinearModel< casacore::Lattice<casacore::Float> >
149 : {
150 :
151 : // Any new entropies derived from Entropy must sign the friend list:
152 : friend class IncEntropy;
153 : friend class IncEntropyI;
154 : friend class IncEntropyEmptiness;
155 :
156 :
157 : public:
158 :
159 : // Construct the CEMemModel object and initialise the model.
160 : IncCEMemModel(IncEntropy &ent,
161 : casacore::Lattice<casacore::Float> & model, // previous model, logically const
162 : casacore::Lattice<casacore::Float> & deltaModel, // this one changes
163 : casacore::uInt nIntegrations = 10,
164 : casacore::Float sigma = 0.001,
165 : casacore::Float targetFlux = 1.0,
166 : casacore::Bool useFluxConstraint = false,
167 : casacore::Bool initializeModel = true,
168 : casacore::Bool imagePlane = false);
169 :
170 : // Construct the CEMemModel object, initialise the model and Prior images.
171 : IncCEMemModel(IncEntropy &ent,
172 : casacore::Lattice<casacore::Float> & model,
173 : casacore::Lattice<casacore::Float> & deltaModel,
174 : casacore::Lattice<casacore::Float> & prior,
175 : casacore::uInt nIntegrations = 10,
176 : casacore::Float sigma = 0.001,
177 : casacore::Float targetFlux = 1.0,
178 : casacore::Bool useFluxConstraint = false,
179 : casacore::Bool initializeModel = true,
180 : casacore::Bool imagePlane = false);
181 :
182 : // Construct the CEMemModel object, initialise the model and
183 : // and mask images.
184 : IncCEMemModel(IncEntropy &ent,
185 : casacore::Lattice<casacore::Float> & model,
186 : casacore::Lattice<casacore::Float> & deltaModel,
187 : casacore::uInt nIntegrations,
188 : casacore::Lattice<casacore::Float> & mask,
189 : casacore::Float sigma = 0.001,
190 : casacore::Float targetFlux = 1.0,
191 : casacore::Bool useFluxConstraint = false,
192 : casacore::Bool initializeModel = true,
193 : casacore::Bool imagePlane = false);
194 :
195 : // Construct the CEMemModel object, initialise the model, Prior,
196 : // and mask images.
197 : IncCEMemModel(IncEntropy &ent,
198 : casacore::Lattice<casacore::Float> & model,
199 : casacore::Lattice<casacore::Float> & deltaModel,
200 : casacore::Lattice<casacore::Float> & prior,
201 : casacore::Lattice<casacore::Float> & mask,
202 : casacore::uInt nIntegrations = 10,
203 : casacore::Float sigma = 0.001,
204 : casacore::Float targetFlux = 1.0,
205 : casacore::Bool useFluxConstraint = false,
206 : casacore::Bool initializeModel = true,
207 : casacore::Bool imagePlane = false);
208 :
209 :
210 : // destructor
211 : virtual ~IncCEMemModel();
212 :
213 :
214 : // solve the convolution equation
215 : // returns true if converged
216 :
217 : // Gives information about the state of the CEMem
218 : // (ie, using mask image, using prior image; more work here!)
219 : void state();
220 :
221 :
222 : // This needs to be "ResidualEquation", using LatConvEquation as
223 : // polymorphism is broken
224 : casacore::Bool solve(ResidualEquation<casacore::Lattice<casacore::Float> > & eqn);
225 :
226 : // Set and get various state images and classes
227 : // <group>
228 : // set or get the Entropy class
229 : void setEntropy(IncEntropy &myEntropy ) {itsEntropy_ptr = &myEntropy;}
230 : void getEntropy(IncEntropy &myEntropy ) {myEntropy = *itsEntropy_ptr;}
231 :
232 : // set or get the Model image
233 0 : casacore::Lattice<casacore::Float>& getModel() const
234 0 : { return (*(itsDeltaModel_ptr)); }
235 0 : void setModel(const casacore::Lattice<casacore::Float> & model)
236 0 : { itsDeltaModel_ptr = model.clone(); }
237 :
238 : // set or get the Prior image
239 : casacore::Lattice<casacore::Float>& getPrior() const
240 : { return (*(itsPrior_ptr->clone())); }
241 : void setPrior(casacore::Lattice<casacore::Float> & prior);
242 :
243 : // set or get the Mask image
244 : casacore::Lattice<casacore::Float>& getMask() const
245 : { return (*(itsMask_ptr->clone())); }
246 : void setMask(casacore::Lattice<casacore::Float> & mask);
247 : // </group>
248 :
249 :
250 : // set and get alpha, beta, and Q
251 : // <group>
252 0 : casacore::Float getAlpha() const { return itsAlpha; }
253 0 : casacore::Float getBeta() const { return itsBeta; }
254 0 : void setAlpha(casacore::Float alpha) {itsAlpha = alpha; }
255 0 : void setBeta(casacore::Float beta) {itsBeta = beta; }
256 : // </group>
257 :
258 : // Set various controlling parameters
259 : // (The most popular controlling parameters are
260 : // set in the constructor)
261 : // <group>
262 : casacore::Float getTolerance() {return itsTolerance; }
263 : void setTolerance(casacore::Float x) { itsTolerance = x; }
264 0 : casacore::Float getQ() {return itsQ; }
265 0 : void setQ(casacore::Float x) { itsQ= x; }
266 : casacore::Float getGain() {return itsGain; }
267 0 : void setGain(casacore::Float x) { itsGain = x; }
268 : casacore::Float getMaxNormGrad() {return itsMaxNormGrad; }
269 : void setMaxNormGrad(casacore::Float x) { itsMaxNormGrad = x; }
270 : casacore::Int getInitialNumberIterations() {return itsFirstIteration; }
271 0 : void setInitialNumberIterations(casacore::Int x) { itsFirstIteration = x; }
272 : // </group>
273 :
274 : // The convergence can also be in terms of the maximum residual
275 : // (ie, for use in stopping in a major cycle).
276 0 : void setThreshold (const casacore::Float x ) { itsThreshold0 = x; }
277 : // Thresh doubles in iter iterations
278 0 : void setThresholdSpeedup (const casacore::Float iter) {itsThresholdSpeedup = iter; }
279 : casacore::Float getThreshold();
280 :
281 : // Set/get the progress display
282 : // <group>
283 0 : virtual void setProgress(CEMemProgress& ccp) { itsProgressPtr = &ccp; }
284 0 : virtual CEMemProgress& getProgress() { return *itsProgressPtr; }
285 : // </group>
286 :
287 : // return the number of iterations completed
288 0 : casacore::Int numberIterations () { return itsIteration; }
289 :
290 : protected:
291 :
292 :
293 : // Perform One Iteration
294 : void oneIteration();
295 :
296 : // apply mask to a lattice; returns true if mask is available,
297 : // false if not
298 : casacore::Bool applyMask( casacore::Lattice<casacore::Float> & lat );
299 :
300 : // Helper functions that interface with Entropy routines
301 : // Access to the entropy should be through this interface; the
302 : // points at which the Entropy is mentioned is then limited to
303 : // these lines right here, and to the constructors which set the
304 : // Entropy. The entropy should not ever change the private
305 : //
306 : // <group>
307 0 : void formEntropy() { itsEntropy = itsEntropy_ptr->formEntropy(); }
308 :
309 0 : void formGDG() { itsEntropy_ptr->formGDG( itsGDG ); }
310 :
311 0 : void formGDGStep() { itsEntropy_ptr->formGDGStep( itsGDG ); }
312 :
313 0 : void formGDS() { itsGradDotStep1=itsEntropy_ptr->formGDS(); }
314 :
315 : void entropyType(casacore::String &str) { itsEntropy_ptr->entropyType(str); }
316 :
317 0 : void relaxMin() { itsRequiredModelMin = itsEntropy_ptr->relaxMin(); }
318 :
319 : casacore::Bool testConvergence() { return itsEntropy_ptr->testConvergence(); }
320 :
321 : // </group>
322 :
323 :
324 : // protected generic constrcutor: DON'T USE IT!
325 : IncCEMemModel();
326 :
327 : // Set entropy pointer to zero: called by Entropy's destructor
328 : void letEntropyDie() { itsEntropy_ptr = 0; }
329 :
330 : // initialize itsStep and itsResidual and other stuff
331 : casacore::Bool initStuff();
332 :
333 :
334 : // controls how to change Alpha and Beta
335 : // <group>
336 : void changeAlphaBeta ();
337 : // initialize Alpha-Beta (called by changeAlphaBeta)
338 : void initializeAlphaBeta();
339 : // update Alpha-Beta (called by changeAlphaBeta)
340 : void updateAlphaBeta();
341 : // </group>
342 :
343 :
344 :
345 :
346 : // Generic utility functions
347 : // <group>
348 : // check that a single image is onf plausible shape
349 : casacore::Bool checkImage(const casacore::Lattice<casacore::Float> *);
350 : // check that the lattices and the underlying tile sizes are consistent
351 : casacore::Bool checkImages(const casacore::Lattice<casacore::Float> *one, const casacore::Lattice<casacore::Float> *other);
352 : // check that all is well in Denmark:
353 : // ensure all images are the same size,
354 : // ensure we have an entropy,
355 : // ensure state variables have reasonable values
356 : casacore::Bool ok();
357 : // </group>
358 :
359 :
360 :
361 :
362 : // Helper functions for oneIteration:
363 : // <group>
364 : // calculate size of step
365 : void calculateStep();
366 :
367 : // take one step: clipped addition of
368 : // wt1*itsModel + wt2*itsStep
369 : void takeStep(casacore::Float wt1, casacore::Float wt2);
370 :
371 : // Calculate the flux, itsModMin, and itsModMax
372 : casacore::Float formFlux();
373 : // </group>
374 :
375 : // Determine if the peak residual is less than the getThreshold()
376 : casacore::Bool testConvergenceThreshold();
377 :
378 : // ------------Private Member casacore::Data---------------------
379 : // functional form of the entropy
380 : IncEntropy *itsEntropy_ptr;
381 :
382 : // form of the Residual method
383 : ResidualEquation< casacore::Lattice<casacore::Float> > * itsResidualEquation_ptr;
384 :
385 : // Images:
386 : casacore::Lattice<casacore::Float> * itsModel_ptr;
387 : casacore::Lattice<casacore::Float> * itsDeltaModel_ptr;
388 : casacore::Lattice<casacore::Float> * itsPrior_ptr;
389 : casacore::Lattice<casacore::Float> * itsMask_ptr;
390 : // Our OWN internal temp images; delete these upon destruction
391 : casacore::Lattice<casacore::Float> * itsStep_ptr;
392 : casacore::Lattice<casacore::Float> * itsResidual_ptr;
393 :
394 :
395 : // Controlling parameters
396 : // <group>
397 : casacore::Bool itsInitializeModel;
398 : casacore::uInt itsNumberIterations;
399 : casacore::Bool itsDoInit;
400 : casacore::Float itsSigma;
401 : casacore::Float itsTargetFlux;
402 : casacore::Float itsDefaultLevel;
403 : casacore::Float itsTargetChisq;
404 : // tolerance for convergence
405 : casacore::Float itsTolerance;
406 : // N points per beam
407 : casacore::Float itsQ;
408 : // gain for adding step image
409 : casacore::Float itsGain;
410 : casacore::Float itsMaxNormGrad;
411 : // constrain flux or not?
412 : casacore::Bool itsUseFluxConstraint;
413 : // is this an image plane problem (like Single Dish or Optical?)
414 : casacore::Bool itsDoImagePlane;
415 : casacore::Float itsThreshold0;
416 : casacore::Float itsThresholdSpeedup;
417 : // </group>
418 :
419 : // State variables
420 : // <group>
421 : casacore::Float itsAlpha;
422 : casacore::Float itsBeta;
423 : casacore::Float itsNormGrad;
424 : // note that itsModelFlux + itsDeltaFlux = itsFlux
425 : casacore::Float itsFlux; // this is the total Flux
426 : casacore::Float itsModelFlux; // this is the flux of image itsModel_ptr
427 : casacore::Float itsDeltaFlux; // this is the flux of image itsDeltaModel_ptr
428 : casacore::Float itsChisq;
429 : // sqrt( chi^2/target_chi^2 )
430 : casacore::Float itsFit;
431 : // sqrt( chi^2/ Npixels )
432 : casacore::Float itsAFit;
433 : // numerical value of entropy
434 : casacore::Float itsEntropy;
435 : // Model is constrained to be >= itsNewMin
436 : casacore::Float itsRequiredModelMin;
437 : // maximum pixel value in model
438 : casacore::Float itsModelMax;
439 : // minimum pixel value n model
440 : casacore::Float itsModelMin;
441 : casacore::Float itsLength;
442 : casacore::Double itsGradDotStep0;
443 : casacore::Double itsGradDotStep1;
444 : casacore::uInt itsIteration;
445 : casacore::uInt itsFirstIteration;
446 : // matrix of gradient dot products
447 : casacore::Matrix<casacore::Double> itsGDG;
448 : casacore::Float itsCurrentPeakResidual;
449 : // </group>
450 : casacore::Float itsNumberPixels;
451 :
452 :
453 : // Accesories
454 : // <group>
455 : casacore::Bool itsChoose;
456 : casacore::LogIO itsLog;
457 : // </group>
458 :
459 : // Enumerate the different Gradient subscript types
460 : enum GradientType {
461 : // Entropy
462 : H=0,
463 : // Chi_sq
464 : C=1,
465 : // Flux
466 : F=2,
467 : // Objective function J
468 : J=3
469 : };
470 :
471 : CEMemProgress* itsProgressPtr;
472 :
473 :
474 : };
475 :
476 :
477 :
478 : } //# NAMESPACE CASA - END
479 :
480 : #endif
481 :
|
