Line data    Source code

       1             : //# ImageSkyModel.h: Definition for ImageSkyModel
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2002
       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 adressed 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$
      28             : 
      29             : #ifndef SYNTHESIS_IMAGESKYMODEL_H
      30             : #define SYNTHESIS_IMAGESKYMODEL_H
      31             : 
      32             : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
      33             : #include <synthesis/MeasurementComponents/SkyModel.h>
      34             : #include <casacore/casa/Arrays/Matrix.h>
      35             : #include <casacore/images/Images/ImageInterface.h>
      36             : #include <casacore/images/Images/PagedImage.h>
      37             : #include <casacore/images/Images/TempImage.h>
      38             : #include <casacore/casa/Logging/LogMessage.h>
      39             : #include <casacore/casa/Logging/LogSink.h>
      40             : #include <casacore/casa/System/PGPlotter.h>
      41             : 
      42             : namespace casa { //# NAMESPACE CASA - BEGIN
      43             : 
      44             : // <summary> 
      45             : // Image Sky Model: Image-based Model for the Sky Brightness
      46             : // </summary>
      47             : 
      48             : // <use visibility=export>
      49             : 
      50             : // <reviewed reviewer="" date="" tests="" demos="">
      51             : 
      52             : // <prerequisite>
      53             : //   <li> <linkto class=SkyModel>SkyModel</linkto> class
      54             : //   <li> <linkto class=SkyEquation>SkyEquation</linkto> class
      55             : //   <li> <linkto class=casacore::ImageInterface>ImageInterface</linkto> class
      56             : //   <li> <linkto class=casacore::PagedImage>PagedImage</linkto> class
      57             : //   <li> <linkto module=MeasurementComponents>MeasurementComponents</linkto> module
      58             : //   <li> <linkto class=VisSet>VisSet</linkto> class
      59             : // </prerequisite>
      60             : //
      61             : // <etymology>
      62             : // ImageSkyModel describes an interface for Models to be used in
      63             : // the SkyEquation. It is derived from <linkto class=SkyModel>SkyModel</linkto>.
      64             : // </etymology>
      65             : //
      66             : // <synopsis> 
      67             : // A ImageSkyModel contains a number of separate models. The interface to
      68             : // SkyEquation is via an image per model. <linkto class=SkyEquation>SkyEquation</linkto> uses this image to
      69             : // calculate Fourier transforms, etc. Some (most) SkyModels are
      70             : // solvable: the SkyEquation can be used by the SkyModel to return
      71             : // gradients with respect to itself (via the image interface). Thus
      72             : // for a SkyModel to solve for itself, it calls the SkyEquation
      73             : // methods to get gradients of chi-squared with respect to the
      74             : // image pixel values (thus returning an image: basically a residual
      75             : // image). The SkyModel then uses these gradients as appropriate to
      76             : // update itself.
      77             : // </synopsis> 
      78             : //
      79             : // <example>
      80             : // See the example for <linkto class=SkyModel>SkyModel</linkto>.
      81             : // </example>
      82             : //
      83             : // <motivation>
      84             : // The properties of a model of the sky must be described
      85             : // for the <linkto class=SkyEquation>SkyEquation</linkto>.
      86             : // </motivation>
      87             : //
      88             : // <todo asof="97/10/01">
      89             : // <li> Multiple images in SkyModel
      90             : // <li> ComponentModel
      91             : // </todo>
      92             : 
      93             : class ImageSkyModel : public SkyModel {
      94             : public:
      95             : 
      96             :   // Empty constructor
      97             :   ImageSkyModel(const casacore::Int maxNumModels=1);
      98             : 
      99             :   void setMaxNumberModels(const casacore::Int maxNumModels);
     100             : 
     101             :   // Copy constructor
     102             :   ImageSkyModel(const ImageSkyModel& sm);
     103             : 
     104             :   // Add a componentlist
     105             :   virtual casacore::Bool add(ComponentList& compList);
     106             :   //update componentlist
     107             :   virtual casacore::Bool updatemodel(ComponentList& compList);
     108             : 
     109             : 
     110             :   // Add an image. maxNumXfr is the maximum Number of transfer functions
     111             :   // that we might want to associate with this image.
     112             :   virtual casacore::Int add(casacore::ImageInterface<casacore::Float>& image, const casacore::Int maxNumXfr=100);
     113             :   //update model image...you have to have added it before...nmodels_p held has to be bigger that image here
     114             :   //its left to the caller to make sure the image is conformant...otherwise you are in trouble.
     115             :   virtual casacore::Bool  updatemodel(const casacore::Int thismodel, casacore::ImageInterface<casacore::Float>& image);
     116             :   // Add a residual image
     117             :   virtual casacore::Bool addResidual(casacore::Int image, casacore::ImageInterface<casacore::Float>& residual);
     118             : 
     119             :   // Destructor
     120             :   virtual ~ImageSkyModel();
     121             : 
     122             :   // Assignment operator
     123             :   ImageSkyModel& operator=(const ImageSkyModel& other);
     124             : 
     125             :   // Number of models contained
     126        1705 :   virtual casacore::Int numberOfModels() {return nmodels_p;};
     127             : 
     128             :   // MFS : Number of taylor terms per model
     129         130 :   virtual casacore::Int numberOfTaylorTerms() {return 1;};
     130             : 
     131             :   // MFS : In-place coefficient residual calculations
     132           0 :   virtual casacore::Bool calculateCoeffResiduals(){return false;};
     133             : 
     134             :   // MFS : Calculate restored alpha and beta.
     135           0 :   virtual casacore::Bool calculateAlphaBeta(const casacore::Vector<casacore::String>& /*restoredNames*/, const casacore::Vector<casacore::String>& /*residualNames*/){return false;};
     136             : 
     137             :   // MFS : Reference Frequency
     138           0 :   virtual casacore::Double getReferenceFrequency(){return 0.0;}
     139             : 
     140             :   // MFS : Index of Taylor term in array of nmodels x ntaylorterms
     141             :   //virtual casacore::Int getTaylorIndex(casacore::Int index){return 0;}
     142           0 :   virtual casacore::Int getTaylorIndex(casacore::Int index){return (casacore::Int)(index/nfields_p);}
     143             : 
     144             :   // Is this model solveable?
     145             :   casacore::Bool isSolveable(casacore::Int model=0);
     146             : 
     147             :   // Free and fix the model (returns previous status). Free means that
     148             :   // it will be solved for in any solution.
     149             :   casacore::Bool free(casacore::Int model=0);
     150             :   casacore::Bool  fix(casacore::Int model=0);
     151             : 
     152             :   // Initialize for gradient search
     153             :   virtual void initializeGradients();
     154             : 
     155             :   // Finalize for gradient search
     156          14 :   virtual void finalizeGradients() {};
     157             : 
     158             :   // Does this have a component list?
     159             :   casacore::Bool hasComponentList();
     160             : 
     161             :   casacore::Bool isEmpty(casacore::Int model=0);
     162             : 
     163             :   // Return the component list
     164             :   virtual ComponentList& componentList();
     165             : 
     166             :   // Return actual images to be used by SkyEquation. 
     167             :   // <group>
     168             :   casacore::ImageInterface<casacore::Float>& image(casacore::Int model=0);
     169             :   casacore::ImageInterface<casacore::Complex>& cImage(casacore::Int model=0);
     170             :   casacore::ImageInterface<casacore::Complex>& XFR(casacore::Int model=0, casacore::Int numXFR=0);
     171             :   casacore::ImageInterface<casacore::Float>& PSF(casacore::Int model=0);
     172             :   casacore::ImageInterface<casacore::Float>& gS(casacore::Int model=0);
     173             :   casacore::ImageInterface<casacore::Float>& residual(casacore::Int model=0);
     174             :   casacore::ImageInterface<casacore::Float>& ggS(casacore::Int model=0);
     175             :   // if (doFluxScale(mod))  image(mod) *  fluxScale(mod)  
     176             :   // gives actual brightness distribution
     177             :   casacore::ImageInterface<casacore::Float>& fluxScale(casacore::Int model=0);
     178             :   casacore::ImageInterface<casacore::Float>& work(casacore::Int model=0);
     179             :   casacore::ImageInterface<casacore::Float>& deltaImage(casacore::Int model=0);
     180             :   // </group>
     181             : 
     182             :   // tells if this model needs to be multiplied by a flux scale image
     183             :   casacore::Bool doFluxScale(casacore::Int model=0);
     184             :   // require use of flux scale image
     185             :   void mandateFluxScale(casacore::Int model=0);
     186             : 
     187             :   casacore::Bool hasXFR(casacore::Int model=0);
     188             : 
     189             :   // Add to Sum weights, Chi-Squared
     190          14 :   void addStatistics(casacore::Float sumwt, casacore::Float chisq) {sumwt_p+=sumwt;chisq_p+=chisq;}
     191             : 
     192             :   // Weight per model (channels, polarizations)
     193             :   casacore::Matrix<casacore::Float>& weight(casacore::Int model=0);
     194             : 
     195             :   // Solve for this SkyModel: This replaces the image with
     196             :   // the residual image
     197             :   virtual casacore::Bool solve (SkyEquation& me);
     198             : 
     199             :   // Solve explicitly for the residuals: same as solve for
     200             :   // this class
     201             :   // modelToMs determines if predicted vis is put in the MODEL_DATA column
     202             :   casacore::Bool solveResiduals (SkyEquation& me, casacore::Bool modelToMS=false);
     203             : 
     204             :   // Make the approximate PSFs needed for each model
     205             :   virtual void makeApproxPSFs(SkyEquation& se);
     206             : 
     207             :   // Get current residual image: this is either that image specified via
     208             :   // addResidual, or a scratch image.
     209             :   // For example in WFImageSkyModel it might return the whole main image
     210             :   //rather than facets 
     211             :   virtual casacore::ImageInterface<casacore::Float>& getResidual (casacore::Int model=0);
     212             : 
     213             :   // Return the fitted beam for each model
     214             :   casacore::ImageBeamSet& beam(casacore::Int model=0);
     215             : 
     216             :   // Set casacore::PGPlotter to be used
     217             :   void setPGPlotter(casacore::PGPlotter& pgp) { pgplotter_p = &pgp; }
     218             : 
     219             :   // This is the factor by which you multiply the worst outer
     220             :   // sidelobe by to get the threshold for the current cycle
     221           0 :   void setCycleFactor(float x) { cycleFactor_p = x; }
     222             : 
     223             :   // Cycle threshold will double in this number of iterations
     224             :   // (ie, use a large number if you don't want cycle threshold
     225             :   // to inch up)
     226           0 :   void setCycleSpeedup(float x) { cycleSpeedup_p = x; }
     227             : 
     228             :   // Yet another control for the minor cycle threshold.
     229             :   // This is in response to CAS-2673
     230             :   // This allows control similar to 'cyclefactor' - used in MFClarkCleanSkyModel
     231           0 :   void setCycleMaxPsfFraction(float x) { cycleMaxPsfFraction_p = x; }
     232             : 
     233             :   // Set the variable that switches on the progress display
     234           0 :   void setDisplayProgress (const casacore::Bool display ) {displayProgress_p = display; };
     235             : 
     236             :   // Set a variable to indicate the polarization frame in the data (circular or linear).
     237             :   // This is used along with the user's choice of output casacore::Stokes parameter
     238             :   // to decide the stokesCoordinate of the temporary images "cImage".
     239         124 :   void setDataPolFrame(StokesImageUtil::PolRep datapolrep) {dataPolRep_p = datapolrep;};
     240             : 
     241             :   // Tries to return a pointer to a casacore::TempImage (allocated with new, so remember
     242             :   // to use delete) with the given shape and CoordinateSystem.
     243             :   //
     244             :   // @param imgShp
     245             :   // @param imgCoords
     246             :   // @param nMouthsToFeed: If > 1 it is taken as a hint that it should leave
     247             :   //                       room for nMouthsToFeed - 1 more TempImages. 
     248             :   //
     249             :   // <throws>
     250             :   // casacore::AipsError on memory allocation error.
     251             :   // </throws>
     252             :   template<class M>
     253             :   static casacore::TempImage<M>* getTempImage(const casacore::TiledShape& imgShp,
     254             :                                     const casacore::CoordinateSystem& imgCoords,
     255             :                                     const casacore::uInt nMouthsToFeed=1);
     256             :   
     257           0 :   virtual casacore::Int getModelIndex(casacore::uInt field, casacore::uInt /*taylor*/){return field;};
     258             : 
     259             :   //try to make templattices use memory if possible
     260             :   //if set to false then always use disk
     261             :   virtual void setMemoryUse(casacore::Bool useMem=false);
     262          14 :   virtual casacore::Bool getMemoryUse(){return useMem_p;};
     263             :   //Set templattice tile vol  in pixels
     264             :   void setTileVol(const casacore::Int tileVol=1000000);
     265             : protected:
     266             : 
     267             :   // Make Newton Raphson step internally. This is really an implementation
     268             :   // detail: it is useful for derived classes.
     269             :   // The modelToMS parameter is for committing to MODEL_DATA column of the MS
     270             :   // the predicted visibilities.
     271             : 
     272             :   casacore::Bool makeNewtonRaphsonStep(SkyEquation& se, 
     273             :                              casacore::Bool incremental=false, casacore::Bool modelToMS=false);
     274             : 
     275             : 
     276             :   // Get casacore::PGPlotter to be used
     277             :   casacore::PGPlotter* getPGPlotter() { return pgplotter_p; }
     278             : 
     279             :   casacore::Int maxnmodels_p;
     280             :   casacore::Int nmodels_p;
     281             :   //MFS
     282             :   casacore::Int nfields_p;
     283             : 
     284             :   casacore::Int maxNumXFR_p;
     285             : 
     286             :   casacore::Float sumwt_p; 
     287             :   casacore::Float chisq_p; 
     288             : 
     289             :   // ComponentList
     290             :   ComponentList* componentList_p;
     291             : 
     292             :   // Images
     293             :   casacore::Vector<casacore::String> imageNames_p;
     294             :   // Everything here can be just interface
     295             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > image_p;
     296             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > residual_p;
     297             : 
     298             :   // We actually create these
     299             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Complex> * > cimage_p;
     300             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Complex> * > cxfr_p;
     301             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > residualImage_p;
     302             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > gS_p;
     303             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > psf_p;
     304             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > ggS_p;
     305             :   // if (doFluxScale_p), image_p * fluxScale_p gives the true brightness
     306             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > fluxScale_p;
     307             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > work_p;
     308             :   casacore::PtrBlock<casacore::ImageInterface<casacore::Float> * > deltaimage_p;
     309             :   casacore::Block<casacore::Bool> solve_p;
     310             :   casacore::Block<casacore::Bool> doFluxScale_p;
     311             : 
     312             :   casacore::PtrBlock<casacore::Matrix<casacore::Float> * > weight_p;
     313             : 
     314             :   casacore::PtrBlock<casacore::ImageBeamSet * > beam_p;
     315             : 
     316             :   casacore::LogSink logSink_p;
     317           0 :   casacore::LogSink& logSink() {return logSink_p;};
     318             :   
     319             :   casacore::Long cacheSize(casacore::Int model);
     320             :   casacore::IPosition tileShape(casacore::Int model);
     321             : 
     322             :   casacore::PGPlotter *pgplotter_p;
     323             :   casacore::Bool displayProgress_p;
     324             :   // This is the factor by which you multiply the worst outer
     325             :   // sidelobe by to get the threshold for the current cycle
     326             :   casacore::Float cycleFactor_p;
     327             :   // Cycle threshold will double in this number of iterations
     328             :   // (ie, use a large number if you don't want cycle threshold
     329             :   // to inch up)
     330             :   casacore::Float cycleSpeedup_p;
     331             :   // Cycle threshold = maxResidual x min(Max-Psf-Fraction , cyclefactor x maxpsfsidelobe)
     332             :   casacore::Float cycleMaxPsfFraction_p;
     333             :   // If PSF is done..should not redo it.
     334             :   casacore::Bool donePSF_p;
     335             :   // check if model has been modified especially for continuing
     336             :   // a deconvolution
     337             :   casacore::Bool modified_p;
     338             :   // Parameter to indicate the polaraization type of the data (circular or linear)
     339             :   // Required by cImage() to decide shapes.
     340             :   StokesImageUtil::PolRep dataPolRep_p;
     341             :   casacore::Bool workDirOnNFS_p;
     342             :   casacore::Bool useMem_p;
     343             :   casacore::Int tileVol_p;
     344             : };
     345             : 
     346             : 
     347             : 
     348             : } //# NAMESPACE CASA - END
     349             : 
     350             : 
     351             : #ifndef AIPS_NO_TEMPLATE_SRC
     352             : #include <synthesis/MeasurementComponents/ImageSkyModel.tcc>
     353             : #endif //# AIPS_NO_TEMPLATE_SRC
     354             : 
     355             : #endif
     356             : 
     357             :