00001 //# ClarkCleanLatModel.h: this defines ClarkCleanLatModel 00002 //# Copyright (C) 1996,1997,1998,1999,2000,2003 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# 00027 //# $Id$ 00028 00029 #ifndef SYNTHESIS_CLARKCLEANLATMODEL_H 00030 #define SYNTHESIS_CLARKCLEANLATMODEL_H 00031 00032 #include <casa/aips.h> 00033 #include <synthesis/MeasurementEquations/LinearModel.h> 00034 #include <synthesis/MeasurementEquations/LinearEquation.h> 00035 #include <lattices/Lattices/Lattice.h> 00036 #include <casa/Arrays/IPosition.h> 00037 #include <synthesis/MeasurementEquations/Iterate.h> 00038 #include <casa/Logging/LogIO.h> 00039 00040 namespace casa { //# NAMESPACE CASA - BEGIN 00041 00042 template <class T> class Matrix; 00043 template <class T> class Vector; 00044 class ClarkCleanProgress; 00045 class LatConvEquation; 00046 class CCList; 00047 00048 // <summary> 00049 // A Class for performing the Clark Clean Algorithm on Arrays 00050 // </summary> 00051 00052 // <use visibility=export> 00053 00054 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos=""> 00055 // </reviewed> 00056 00057 // <prerequisite> 00058 // <li> ResidualEquation/LatConvEquation 00059 // <li> LinearModel/LinearEquation Paradigm 00060 // </prerequisite> 00061 // 00062 // <etymology> 00063 // This class is called ClarkCleanLatModel because thats the algorithm it uses 00064 // deconvolve the lattice-based model. 00065 // </etymology> 00066 // 00067 // <synopsis> 00068 // This class is used to perform the Clark Clean Algorithm on an 00069 // Array. Only the deconvolved model of the sky are directly stored by this 00070 // class. The point spread function (psf) and convolved (dirty) image are 00071 // stored in a companion class which is must be derived from 00072 // ResidualEquation. 00073 // 00074 // The cleaning works like this. The user constructs a ClarkCleanLatModel by 00075 // specifying an initial model of the sky. This can by be 00076 // one,two,three... dimensional depending on the dimension of the psf (see 00077 // below). The user then constructs a class which implements the forward 00078 // equation between the model and the dirty image. Typically this will be 00079 // the ConvolutionEquation class, although any class which has a 00080 // ResidualEquation interface will be work (but perhaps very slowly, as the 00081 // ConvolutionEquation class has member functions optimised for cleaning) 00082 // 00083 // The user then calls the solve() function (with the appropriate equation 00084 // class as an arguement), and this class will perform the Clark clean. 00085 // The various clean parameters are set (prior to calling solve) using the 00086 // functions derived from the Iterate class, in particular setGain(), 00087 // setNumberIterations() & setThreshold() (to set a flux limit). 00088 // 00089 // The solve() function does not return either the deconvolved model or the 00090 // residuals. The solved model can be obtained using the getModel() function 00091 // (derived from ArrayModel()) and the residual can be obtained using the 00092 // residual() member function of the Convolution/Residual Equation Class. 00093 // 00094 // The size and shape of the model used in this class MUST be the same as 00095 // the convolved data (Dirty Image), stored in the companion 00096 // ResidualEquation Class. However the model (and convolved data) can have 00097 // more dimensions than the psf, as well as a different size (either larger 00098 // or smaller). When the dimensionality is different the cleaning is done 00099 // independendtly in each "plane" of the model. (Note this has not 00100 // been implemented yet but is relatively simple to do if necessary). 00101 // 00102 // This multi-dimensionalty is exploited when cleaning arrays of 00103 // StokesVectors. Here the Array of StokesVectors is decomposed into a stack 00104 // of 4 Floating point arrays and the cleaning is done on all the the arrays 00105 // simultaneosly. The criterion for choosing the brightest pixel has been 00106 // generalised by using the "length" of the Stokesvector in 4 dimensional 00107 // space. 00108 // 00109 // A companion class to this one is MaskedClarkCleanLatModel. This provides 00110 // the same functionality but is used with MaskedArrays which indicate which 00111 // regions of the model to search for clean components. 00112 // 00113 // </synopsis> 00114 // 00115 // <example> 00116 // <srcblock> 00117 // Matrix<Float> psf(12,12), dirty(10,10), initialModel(10,10); 00118 // ...put appropriate values into psf, dirty, & initialModel.... 00119 // ClarkCleanLatModel<Float> deconvolvedModel(initialModel); 00120 // ConvolutionEquation convEqn(psf, dirty); 00121 // deconvolvedModel.setGain(0.2); 00122 // deconvolvedModel.setNumberIterations(1000); 00123 // Bool convWorked = deconvolvedModel.solve(convEqn); 00124 // if (convWorked) 00125 // ConvEqn.residual(deconvolvedModel, finalResidual); 00126 // </srcblock> 00127 // </example> 00128 // 00129 // <motivation> 00130 // This class is needed to deconvolve images. 00131 // </motivation> 00132 // 00133 // <templating arg=T> 00134 // I have tested this class with Arrays of 00135 // <li> Float 00136 // <li> StokesVector 00137 // </templating> 00138 // 00139 // <todo asof="1996/05/02"> 00140 // <li> Make changes so that multidimensions work as advertised 00141 // <li> compare timing with other clean implementations (ie, Mark's 00142 // CleanTools, SDE, AIPS & miriad) 00143 // </todo> 00144 00145 class ClarkCleanLatModel: 00146 public LinearModel< Lattice<Float> >, 00147 public Iterate 00148 { 00149 public: 00150 // The default constructor does nothing more than initialise a zero length 00151 // array to hold the deconvolved model. If this constructor is used then 00152 // the actual model must be set using the setModel() function of the 00153 // LatticeModel class. 00154 ClarkCleanLatModel(); 00155 00156 // Construct the ClarkCleanLatModel object and initialise the model. 00157 ClarkCleanLatModel(Lattice<Float> & model); 00158 00159 // Construct the ClarkCleanLatModel object and initialise the model ans mask 00160 ClarkCleanLatModel(Lattice<Float> & model, Lattice<Float> & mask); 00161 00162 // Construct the ClarkCleanLatModel object and initialise the model ans mask 00163 ClarkCleanLatModel(Lattice<Float> & model, Lattice<Float> & residual, 00164 Lattice<Float> & mask); 00165 // Destroy! 00166 virtual ~ClarkCleanLatModel(); 00167 00168 virtual const Lattice<Float> & getModel() const { return *itsModelPtr; } 00169 virtual void setModel(const Lattice<Float> & model); 00170 virtual void setModel(Lattice<Float> & model); 00171 00172 const Lattice<Float> & getMask() const; 00173 void setMask(const Lattice<Float> & mask); 00174 00175 void setResidual( Lattice<Float> & residual); 00176 virtual const Lattice<Float> & getResidual() const { return *itsResidualPtr; } 00177 00178 Int getNumberIterations(){ return numberIterations(); } 00179 00180 Float getMaxResidual() { return itsMaxRes;}; 00181 // Using a Clark clean deconvolution proceedure solve for an improved 00182 // estimate of the deconvolved object. The convolution/residual equation 00183 // contains the psf and dirty image. When called with a ResidualEquation 00184 // arguement a quite general interface is used that is slow. The 00185 // convolution equation contains functions that speed things up. The 00186 // functions return False if the deconvolution could not be done. 00187 // <group> 00188 Bool solve(LatConvEquation & eqn); 00189 Bool singleSolve(LatConvEquation & eqn, Lattice<Float> & residual); 00190 // </group> 00191 00192 // The user can be asked whether to stop after every minor cycle 00193 // <group> 00194 virtual void setChoose(const Bool askForChoice); 00195 virtual Bool getChoose(); 00196 // </group> 00197 00198 // These remaining functions set various "knobs" that the user can tweak and 00199 // are specific to the Clark clean algorithm. The more generic parameters 00200 // ie. clean gain, and maximum residual fluxlimit, are set using functions in 00201 // the Iterate base class. The get functions return the value that was 00202 // actually used after the cleaning was done. 00203 00204 // set the size of the PSF used in the minor iterations. If not set it 00205 // defaults to the largest useful Psf (ie. min(modelSize*2, psfSize)) 00206 // <group> 00207 virtual void setPsfPatchSize(const IPosition & psfPatchSize); 00208 virtual IPosition getPsfPatchSize(); 00209 // </group> 00210 00211 // Set the size of the histogram used to determine how many pixels are 00212 // "active" in a minor iteration. Default value of 1000 is OK for 00213 // everything except very small cleans. 00214 // <group> 00215 virtual void setHistLength(const uInt histBins); 00216 virtual uInt getHistLength(); 00217 // </group> 00218 00219 // Set the maximum number of minor iterations to perform for each major 00220 // cycle. 00221 // <group> 00222 virtual void setMaxNumberMinorIterations(const uInt maxNumMinorIterations); 00223 virtual uInt getMaxNumberMinorIterations(); 00224 // </group> 00225 00226 // Set and get the initial number of iterations 00227 // <group> 00228 virtual void setInitialNumberIterations(const uInt initialNumberIterations); 00229 virtual uInt getInitialNumberIterations(); 00230 // </group> 00231 00232 // Set the maximum number of major cycles to perform 00233 // <group> 00234 virtual void setMaxNumberMajorCycles(const uInt maxNumMajorCycles); 00235 virtual uInt getMaxNumberMajorCycles(); 00236 // </group> 00237 00238 // Set the maximum number of active pixels to use in the minor 00239 // iterations. The specified number can be exceeded if the topmost bin of 00240 // the histogram contains more pixels than specified here. The default is 00241 // 10,000 which is suitable for images of 512by512 pixels. Reduce this for 00242 // smaller images and increase it for larger ones. 00243 // <group> 00244 virtual void setMaxNumPix(const uInt maxNumPix ); 00245 virtual uInt getMaxNumPix(); 00246 // </group> 00247 00248 00249 // Set the maximum exterior psf value. This is used to determine when to 00250 // stop the minor itartions. This is normally determined from the psf and 00251 // the number set here is only used if this cannot be determined. The 00252 // default is zero. 00253 // <group> 00254 virtual void setMaxExtPsf(const Float maxExtPsf ); 00255 virtual Float getMaxExtPsf(); 00256 // </group> 00257 00258 // The total flux density in the model. 00259 Float modelFlux(); 00260 00261 // An exponent on the F(m,n) factor (see Clark[1980]) which influences how 00262 // quickly active pixels are treated as unreliable. Larger values mean 00263 // more major iterations. The default is zero. I have no experience on 00264 // when to use this factor. 00265 // <group> 00266 virtual void setSpeedup(const Float speedup); 00267 virtual Float getSpeedup(); 00268 // </group> 00269 //Set the cycle factor....the larger this is the shallower is the minor 00270 //cycle 00271 virtual void setCycleFactor(const Float factor); 00272 00273 00274 // Set/get the progress display 00275 // <group> 00276 virtual void setProgress(ClarkCleanProgress& ccp) { itsProgressPtr = &ccp; } 00277 virtual ClarkCleanProgress& getProgress() { return *itsProgressPtr; } 00278 // </group> 00279 00280 private: 00281 // Do all the minor iterations for one major cycle. Cleaning stops 00282 // when the flux or iteration limit is reached. 00283 void doMinorIterations(CCList & activePixels, 00284 Matrix<Float> & psfPatch, 00285 Float fluxLimit, 00286 uInt & numberIterations, 00287 Float Fmn, 00288 const uInt totalIterations, 00289 Float& totalFlux); 00290 00291 void cacheActivePixels(CCList & activePixels, 00292 const Lattice<Float> & residual, Float fluxLimit); 00293 00294 // make histogram of absolute values in array 00295 void absHistogram(Vector<Int> & hist, Float & minVal, 00296 Float & maxVal, const Lattice<Float> & data); 00297 00298 // Determine the flux limit if we only select the maxNumPix biggest 00299 // residuals. Flux limit is not exact due to quantising by the histogram 00300 Float biggestResiduals(Float & maxRes, const uInt maxNumPix, 00301 const Float fluxLimit, 00302 const Lattice<Float> & residual); 00303 00304 // Work out the size of the Psf patch to use. 00305 Float getPsfPatch(Matrix<Float> & psfPatch, LatConvEquation & eqn); 00306 00307 // The maximum residual is the absolute maximum. 00308 Float maxResidual(const Lattice<Float> & residual); 00309 void maxVect(Block<Float> & maxVal, Float & absVal, Int & offset, 00310 const CCList & activePixels); 00311 void subtractComponent(CCList & activePixels, const Block<Float> & maxVal, 00312 const Block<Int> & maxPos, const Matrix<Float> & psf); 00313 Float absMaxBeyondDist(const IPosition & maxDist, const IPosition & centre, 00314 const Lattice<Float> & psf); 00315 Bool stopnow(); 00316 Int getbig(Float const * pixValPtr, Int const * pixPosPtr, const Int nPix, 00317 const Float fluxLimit, 00318 const Float * const residualPtr, const Float * const maskPtr, 00319 const uInt npol, const Int nx, const Int ny); 00320 00321 void updateModel(CCList & cleanComponents); 00322 00323 Lattice<Float> * itsModelPtr; 00324 Lattice<Float> * itsResidualPtr; 00325 const Lattice<Float> * itsSoftMaskPtr; 00326 uInt itsMaxNumPix; 00327 uInt itsHistBins; 00328 uInt itsMaxNumberMinorIterations; 00329 uInt itsInitialNumberIterations; 00330 Int itsMaxNumberMajorCycles; 00331 Float itsMaxExtPsf; 00332 Float itsMaxRes; 00333 IPosition itsPsfPatchSize; 00334 Bool itsChoose; 00335 Float itsSpeedup; 00336 Float itsCycleFactor; 00337 LogIO itsLog; 00338 ClarkCleanProgress* itsProgressPtr; 00339 Bool itsJustStarting; 00340 Bool itsWarnFlag; 00341 }; 00342 00343 00344 } //# NAMESPACE CASA - END 00345 00346 #endif