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