casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
ImageSkyModel.h
Go to the documentation of this file.
00001 //# ImageSkyModel.h: Definition for ImageSkyModel
00002 //# Copyright (C) 1996,1997,1998,1999,2000,2002
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 adressed 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_IMAGESKYMODEL_H
00030 #define SYNTHESIS_IMAGESKYMODEL_H
00031 
00032 #include <ms/MeasurementSets/MeasurementSet.h>
00033 #include <synthesis/MeasurementComponents/SkyModel.h>
00034 #include <casa/Arrays/Matrix.h>
00035 #include <images/Images/ImageInterface.h>
00036 #include <images/Images/PagedImage.h>
00037 #include <images/Images/TempImage.h>
00038 #include <casa/Logging/LogMessage.h>
00039 #include <casa/Logging/LogSink.h>
00040 #include <casa/System/PGPlotter.h>
00041 
00042 namespace casa { //# NAMESPACE CASA - BEGIN
00043 
00044 // <summary> 
00045 // Image Sky Model: Image-based Model for the Sky Brightness
00046 // </summary>
00047 
00048 // <use visibility=export>
00049 
00050 // <reviewed reviewer="" date="" tests="" demos="">
00051 
00052 // <prerequisite>
00053 //   <li> <linkto class=SkyModel>SkyModel</linkto> class
00054 //   <li> <linkto class=SkyEquation>SkyEquation</linkto> class
00055 //   <li> <linkto class=ImageInterface>ImageInterface</linkto> class
00056 //   <li> <linkto class=PagedImage>PagedImage</linkto> class
00057 //   <li> <linkto module=MeasurementComponents>MeasurementComponents</linkto> module
00058 //   <li> <linkto class=VisSet>VisSet</linkto> class
00059 // </prerequisite>
00060 //
00061 // <etymology>
00062 // ImageSkyModel describes an interface for Models to be used in
00063 // the SkyEquation. It is derived from <linkto class=SkyModel>SkyModel</linkto>.
00064 // </etymology>
00065 //
00066 // <synopsis> 
00067 // A ImageSkyModel contains a number of separate models. The interface to
00068 // SkyEquation is via an image per model. <linkto class=SkyEquation>SkyEquation</linkto> uses this image to
00069 // calculate Fourier transforms, etc. Some (most) SkyModels are
00070 // solvable: the SkyEquation can be used by the SkyModel to return
00071 // gradients with respect to itself (via the image interface). Thus
00072 // for a SkyModel to solve for itself, it calls the SkyEquation
00073 // methods to get gradients of chi-squared with respect to the
00074 // image pixel values (thus returning an image: basically a residual
00075 // image). The SkyModel then uses these gradients as appropriate to
00076 // update itself.
00077 // </synopsis> 
00078 //
00079 // <example>
00080 // See the example for <linkto class=SkyModel>SkyModel</linkto>.
00081 // </example>
00082 //
00083 // <motivation>
00084 // The properties of a model of the sky must be described
00085 // for the <linkto class=SkyEquation>SkyEquation</linkto>.
00086 // </motivation>
00087 //
00088 // <todo asof="97/10/01">
00089 // <li> Multiple images in SkyModel
00090 // <li> ComponentModel
00091 // </todo>
00092 
00093 class ImageSkyModel : public SkyModel {
00094 public:
00095 
00096   // Empty constructor
00097   ImageSkyModel(const Int maxNumModels=1);
00098 
00099   void setMaxNumberModels(const Int maxNumModels);
00100 
00101   // Copy constructor
00102   ImageSkyModel(const ImageSkyModel& sm);
00103 
00104   // Add a componentlist
00105   virtual Bool add(ComponentList& compList);
00106   //update componentlist
00107   virtual Bool updatemodel(ComponentList& compList);
00108 
00109 
00110   // Add an image. maxNumXfr is the maximum Number of transfer functions
00111   // that we might want to associate with this image.
00112   virtual Int add(ImageInterface<Float>& image, const Int maxNumXfr=100);
00113   //update model image...you have to have added it before...nmodels_p held has to be bigger that image here
00114   //its left to the caller to make sure the image is conformant...otherwise you are in trouble.
00115   virtual Bool  updatemodel(const Int thismodel, ImageInterface<Float>& image);
00116   // Add a residual image
00117   virtual Bool addResidual(Int image, ImageInterface<Float>& residual);
00118 
00119   // Destructor
00120   virtual ~ImageSkyModel();
00121 
00122   // Assignment operator
00123   ImageSkyModel& operator=(const ImageSkyModel& other);
00124 
00125   // Number of models contained
00126   virtual Int numberOfModels() {return nmodels_p;};
00127 
00128   // MFS : Number of taylor terms per model
00129   virtual Int numberOfTaylorTerms() {return 1;};
00130 
00131   // MFS : In-place coefficient residual calculations
00132   virtual Bool calculateCoeffResiduals(){return False;};
00133 
00134   // MFS : Calculate restored alpha and beta.
00135   virtual Bool calculateAlphaBeta(const Vector<String>& /*restoredNames*/, const Vector<String>& /*residualNames*/){return False;};
00136 
00137   // MFS : Reference Frequency
00138   virtual Double getReferenceFrequency(){return 0.0;}
00139 
00140   // MFS : Index of Taylor term in array of nmodels x ntaylorterms
00141   //virtual Int getTaylorIndex(Int index){return 0;}
00142   virtual Int getTaylorIndex(Int index){return (Int)(index/nfields_p);}
00143 
00144   // Is this model solveable?
00145   Bool isSolveable(Int model=0);
00146 
00147   // Free and fix the model (returns previous status). Free means that
00148   // it will be solved for in any solution.
00149   Bool free(Int model=0);
00150   Bool  fix(Int model=0);
00151 
00152   // Initialize for gradient search
00153   virtual void initializeGradients();
00154 
00155   // Finalize for gradient search
00156   virtual void finalizeGradients() {};
00157 
00158   // Does this have a component list?
00159   Bool hasComponentList();
00160 
00161   Bool isEmpty(Int model=0);
00162 
00163   // Return the component list
00164   virtual ComponentList& componentList();
00165 
00166   // Return actual images to be used by SkyEquation. 
00167   // <group>
00168   ImageInterface<Float>& image(Int model=0);
00169   ImageInterface<Complex>& cImage(Int model=0);
00170   ImageInterface<Complex>& XFR(Int model=0, Int numXFR=0);
00171   ImageInterface<Float>& PSF(Int model=0);
00172   ImageInterface<Float>& gS(Int model=0);
00173   ImageInterface<Float>& residual(Int model=0);
00174   ImageInterface<Float>& ggS(Int model=0);
00175   // if (doFluxScale(mod))  image(mod) *  fluxScale(mod)  
00176   // gives actual brightness distribution
00177   ImageInterface<Float>& fluxScale(Int model=0);
00178   ImageInterface<Float>& work(Int model=0);
00179   ImageInterface<Float>& deltaImage(Int model=0);
00180   // </group>
00181 
00182   // tells if this model needs to be multiplied by a flux scale image
00183   Bool doFluxScale(Int model=0);
00184   // require use of flux scale image
00185   void mandateFluxScale(Int model=0);
00186 
00187   Bool hasXFR(Int model=0);
00188 
00189   // Add to Sum weights, Chi-Squared
00190   void addStatistics(Float sumwt, Float chisq) {sumwt_p+=sumwt;chisq_p+=chisq;}
00191 
00192   // Weight per model (channels, polarizations)
00193   Matrix<Float>& weight(Int model=0);
00194 
00195   // Solve for this SkyModel: This replaces the image with
00196   // the residual image
00197   virtual Bool solve (SkyEquation& me);
00198 
00199   // Solve explicitly for the residuals: same as solve for
00200   // this class
00201   // modelToMs determines if predicted vis is put in the MODEL_DATA column
00202   Bool solveResiduals (SkyEquation& me, Bool modelToMS=False);
00203 
00204   // Make the approximate PSFs needed for each model
00205   virtual void makeApproxPSFs(SkyEquation& se);
00206 
00207   // Get current residual image: this is either that image specified via
00208   // addResidual, or a scratch image.
00209   // For example in WFImageSkyModel it might return the whole main image
00210   //rather than facets 
00211   virtual ImageInterface<Float>& getResidual (Int model=0);
00212 
00213   // Return the fitted beam for each model
00214   ImageBeamSet& beam(Int model=0);
00215 
00216   // Set PGPlotter to be used
00217   void setPGPlotter(PGPlotter& pgp) { pgplotter_p = &pgp; }
00218 
00219   // This is the factor by which you multiply the worst outer
00220   // sidelobe by to get the threshold for the current cycle
00221   void setCycleFactor(float x) { cycleFactor_p = x; }
00222 
00223   // Cycle threshold will double in this number of iterations
00224   // (ie, use a large number if you don't want cycle threshold
00225   // to inch up)
00226   void setCycleSpeedup(float x) { cycleSpeedup_p = x; }
00227 
00228   // Yet another control for the minor cycle threshold.
00229   // This is in response to CAS-2673
00230   // This allows control similar to 'cyclefactor' - used in MFClarkCleanSkyModel
00231   void setCycleMaxPsfFraction(float x) { cycleMaxPsfFraction_p = x; }
00232 
00233   // Set the variable that switches on the progress display
00234   void setDisplayProgress (const Bool display ) {displayProgress_p = display; };
00235 
00236   // Set a variable to indicate the polarization frame in the data (circular or linear).
00237   // This is used along with the user's choice of output Stokes parameter
00238   // to decide the stokesCoordinate of the temporary images "cImage".
00239   void setDataPolFrame(StokesImageUtil::PolRep datapolrep) {dataPolRep_p = datapolrep;};
00240 
00241   // Tries to return a pointer to a TempImage (allocated with new, so remember
00242   // to use delete) with the given shape and CoordinateSystem.
00243   //
00244   // @param imgShp
00245   // @param imgCoords
00246   // @param nMouthsToFeed: If > 1 it is taken as a hint that it should leave
00247   //                       room for nMouthsToFeed - 1 more TempImages. 
00248   //
00249   // <throws>
00250   // AipsError on memory allocation error.
00251   // </throws>
00252   template<class M>
00253   static TempImage<M>* getTempImage(const TiledShape& imgShp,
00254                                     const CoordinateSystem& imgCoords,
00255                                     const uInt nMouthsToFeed=1);
00256   
00257   virtual Int getModelIndex(uInt field, uInt /*taylor*/){return field;};
00258 
00259 protected:
00260 
00261   // Make Newton Raphson step internally. This is really an implementation
00262   // detail: it is useful for derived classes.
00263   // The modelToMS parameter is for committing to MODEL_DATA column of the MS
00264   // the predicted visibilities.
00265 
00266   Bool makeNewtonRaphsonStep(SkyEquation& se, 
00267                              Bool incremental=False, Bool modelToMS=False);
00268 
00269 
00270   // Get PGPlotter to be used
00271   PGPlotter* getPGPlotter() { return pgplotter_p; }
00272 
00273   Int maxnmodels_p;
00274   Int nmodels_p;
00275   //MFS
00276   Int nfields_p;
00277 
00278   Int maxNumXFR_p;
00279 
00280   Float sumwt_p; 
00281   Float chisq_p; 
00282 
00283   // ComponentList
00284   ComponentList* componentList_p;
00285 
00286   // Images
00287   Vector<String> imageNames_p;
00288   // Everything here can be just interface
00289   PtrBlock<ImageInterface<Float> * > image_p;
00290   PtrBlock<ImageInterface<Float> * > residual_p;
00291 
00292   // We actually create these
00293   PtrBlock<ImageInterface<Complex> * > cimage_p;
00294   PtrBlock<ImageInterface<Complex> * > cxfr_p;
00295   PtrBlock<ImageInterface<Float> * > residualImage_p;
00296   PtrBlock<ImageInterface<Float> * > gS_p;
00297   PtrBlock<ImageInterface<Float> * > psf_p;
00298   PtrBlock<ImageInterface<Float> * > ggS_p;
00299   // if (doFluxScale_p), image_p * fluxScale_p gives the true brightness
00300   PtrBlock<ImageInterface<Float> * > fluxScale_p;
00301   PtrBlock<ImageInterface<Float> * > work_p;
00302   PtrBlock<ImageInterface<Float> * > deltaimage_p;
00303   Block<Bool> solve_p;
00304   Block<Bool> doFluxScale_p;
00305 
00306   PtrBlock<Matrix<Float> * > weight_p;
00307 
00308   PtrBlock<ImageBeamSet * > beam_p;
00309 
00310   LogSink logSink_p;
00311   LogSink& logSink() {return logSink_p;};
00312   
00313   Long cacheSize(Int model);
00314 
00315   PGPlotter *pgplotter_p;
00316   Bool displayProgress_p;
00317   // This is the factor by which you multiply the worst outer
00318   // sidelobe by to get the threshold for the current cycle
00319   Float cycleFactor_p;
00320   // Cycle threshold will double in this number of iterations
00321   // (ie, use a large number if you don't want cycle threshold
00322   // to inch up)
00323   Float cycleSpeedup_p;
00324   // Cycle threshold = maxResidual x min(Max-Psf-Fraction , cyclefactor x maxpsfsidelobe)
00325   Float cycleMaxPsfFraction_p;
00326   // If PSF is done..should not redo it.
00327   Bool donePSF_p;
00328   // check if model has been modified especially for continuing
00329   // a deconvolution
00330   Bool modified_p;
00331   // Parameter to indicate the polaraization type of the data (circular or linear)
00332   // Required by cImage() to decide shapes.
00333   StokesImageUtil::PolRep dataPolRep_p;
00334   Bool workDirOnNFS_p;
00335 };
00336 
00337 
00338 
00339 } //# NAMESPACE CASA - END
00340 
00341 
00342 #ifndef AIPS_NO_TEMPLATE_SRC
00343 #include <synthesis/MeasurementComponents/ImageSkyModel.tcc>
00344 #endif //# AIPS_NO_TEMPLATE_SRC
00345 
00346 #endif
00347 
00348