Line data    Source code

       1             : //# WBCleanImageSkyModel.h: Definition for WBCleanImageSkyModel
       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 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_WBCLEANIMAGESKYMODEL_H
      30             : #define SYNTHESIS_WBCLEANIMAGESKYMODEL_H
      31             : 
      32             : #include <synthesis/MeasurementComponents/CleanImageSkyModel.h>
      33             : #include <casacore/lattices/LatticeMath/LatticeCleanProgress.h>
      34             : #include <casacore/lattices/LEL/LatticeExprNode.h>
      35             : #include <casacore/lattices/Lattices/LatticeIterator.h>
      36             : #include <synthesis/MeasurementEquations/MultiTermMatrixCleaner.h>
      37             : #include <casacore/casa/OS/Timer.h>
      38             : 
      39             : namespace casa { //# NAMESPACE CASA - BEGIN
      40             : 
      41             : //forward
      42             : class SkyEquation;
      43             : 
      44             : // <summary> 
      45             : // WB Clean Image Sky Model: Image Sky Model implementing a Wide-Band
      46             : // multi frequency synthesis algorithm
      47             : // </summary>
      48             : 
      49             : // <use visibility=export>
      50             : 
      51             : // <reviewed reviewer="" date="" tests="" demos="">
      52             : 
      53             : // <prerequisite>
      54             : //   <li> <linkto class=casacore::LatticeCleaner>LatticeCleaner</linkto> module
      55             : //   <li> <linkto class=ImageSkyModel>ImageSkyModel</linkto> module
      56             : //   <li> <linkto class=LinearModel>LinearModel</linkto> module
      57             : // </prerequisite>
      58             : //
      59             : // <etymology>
      60             : // WBCleanImageSkyModel implements the Wide Band Clean algorithm.
      61             : // It is derived from <linkto class=SkyModel>SkyModel</linkto>.
      62             : // </etymology>
      63             : //
      64             : // <synopsis> 
      65             : // The WB Clean is the multi-frequency synthesis deconvolution
      66             : // algorithm. It decomposes an image into a linear compbination
      67             : // of models convolved with spectral dirty beams of various
      68             : // order. A multiscale variant can be invoked by supplying a
      69             : // user vector of scale sizes. Default is 1, corresponding to a
      70             : // scale insensitive mfs deconvolution.
      71             : //
      72             : // Masking is optionally performed using a mask image: only points
      73             : // where the mask is non-zero are searched for Gaussian components. 
      74             : // This can cause some difficulty, as the different Gaussian scale
      75             : // sizes will extend beyond the mask by different amounts.
      76             : // If no mask is specified
      77             : // all points in the inner quarter of the image are cleaned.
      78             : // </synopsis> 
      79             : //
      80             : // <example>
      81             : // See the example for <linkto class=SkyModel>SkyModel</linkto>.
      82             : // </example>
      83             : //
      84             : // <motivation>
      85             : // </motivation>
      86             : //
      87             : // <todo asof="99/04/07">
      88             : // <ul> Improve the masking.
      89             : // </todo>
      90             : 
      91             : class WBCleanImageSkyModel : public CleanImageSkyModel {
      92             : public:
      93             : 
      94             :   // Create a WBCleanImageSkyModel - default scale size = 1 pixel
      95             :   WBCleanImageSkyModel();
      96             :   WBCleanImageSkyModel(const casacore::Int ntaylor,const casacore::Int nscales,const casacore::Double reffreq);
      97             :   WBCleanImageSkyModel(const casacore::Int ntaylor,const casacore::Vector<casacore::Float>& userScaleSizes, const casacore::Double reffreq);
      98             : 
      99             :   // destructor
     100             :   ~WBCleanImageSkyModel();
     101             : 
     102             :   // Solve for this SkyModel
     103             :   casacore::Bool solve (SkyEquation& se);
     104             :   //  casacore::Bool copyLatToImInt(casacore::TempLattice<casacore::Float>& lat, casacore::ImageInterface<casacore::Float>& im);
     105             :   // casacore::Bool copyImIntToLat(casacore::TempLattice<casacore::Float>& lat, casacore::ImageInterface<casacore::Float>& im);
     106             :   
     107             : //  casacore::Int nmodels_p; // Number of image models = nfields * ntaylor
     108             :   casacore::Int ntaylor_p; // Number of terms in the Taylor expansion to use.
     109             : //  casacore::Int nfields_p; // Number of image fields/pointings.
     110             :   casacore::Int nscales_p; // Number of scales to use for the multiscale part.
     111             :   
     112             :   casacore::Double refFrequency_p;
     113             : 
     114             : //   casacore::Int add(casacore::ImageInterface<casacore::Float>& iimage, const casacore::Int maxNumXfr=100);
     115             : //   casacore::Bool addResidual(casacore::Int thismodel, casacore::ImageInterface<casacore::Float>& iresidual);
     116             : //   void initializeGradients();
     117             :    casacore::Bool solveResiduals(SkyEquation& se, casacore::Bool modelToMS=false);
     118             :    casacore::Bool makeNewtonRaphsonStep(SkyEquation& se, casacore::Bool incremental=false, casacore::Bool modelToMS=false);
     119             : 
     120           0 :    casacore::Int numberOfTaylorTerms(){return ntaylor_p;};
     121           0 :    casacore::Double getReferenceFrequency(){return refFrequency_p;};
     122             : 
     123             :    // Major axis for ordering : Taylor
     124           0 :    casacore::Int getModelIndex(casacore::uInt model, casacore::uInt taylor){return taylor * (nfields_p) + (model);};
     125             :    casacore::Int getTaylorIndex(casacore::uInt index){return casacore::Int(index/nfields_p);};
     126           0 :    casacore::Int getFieldIndex(casacore::uInt index){return index%nfields_p;};
     127             : 
     128             :   casacore::Bool calculateCoeffResiduals();
     129             :   casacore::Bool calculateAlphaBeta(const casacore::Vector<casacore::String> &restoredNames, 
     130             :                                        const casacore::Vector<casacore::String> &residualNames);
     131             : 
     132             :    // Major axis for ordering : Models
     133             :    //inline casacore::Int getModelIndex(casacore::uInt model, casacore::uInt taylor){return model * (ntaylor_p) + (taylor);};
     134             :    //inline casacore::Int getPSFModelIndex(casacore::uInt model, casacore::uInt taylor){return model * (2*ntaylor_p-1) + (taylor);};
     135             :    //inline casacore::Int getTaylorIndex(casacore::uInt index){return index%ntaylor_p;};
     136             :    //inline casacore::Int getFieldIndex(casacore::uInt index){return index/ntaylor_p;};
     137             :  
     138             :    casacore::Vector<casacore::String> imageNames;
     139             :    
     140             : private:
     141             : 
     142             :   //  casacore::PtrBlock<casacore::MultiTermLatticeCleaner<casacore::Float>* > lc_p;
     143             :   casacore::Block<MultiTermMatrixCleaner> lc_p;
     144             :    
     145             :   casacore::Vector<casacore::Float> scaleSizes_p; // casacore::Vector of scale sizes in pixels.
     146             :   casacore::Vector<casacore::Float> scaleBias_p; // casacore::Vector of scale biases !!
     147             :   //casacore::Float maxPsf_p;
     148             : 
     149             :   casacore::IPosition gip,imshape;
     150             :   casacore::Bool donePSF_p;
     151             :   casacore::Bool doneMTMCinit_p;
     152             :   casacore::Int nx,ny;
     153             : 
     154             :   casacore::Int numbermajorcycles_p;
     155             :   casacore::Float previous_maxresidual_p;
     156             :   
     157             :   // casacore::Memory to be allocated per TempLattice
     158             :   casacore::Double memoryMB_p;
     159             :   
     160             :   casacore::LogIO os;
     161             :   
     162             :   void initVars();
     163             :   casacore::Bool checkParameters();
     164             : 
     165             :   casacore::Int storeAsImg(casacore::String fileName, casacore::ImageInterface<casacore::Float>& theImg);
     166             :   //casacore::Int storeTLAsImg(casacore::String fileName, casacore::TempLattice<casacore::Float> &TL, casacore::ImageInterface<casacore::Float>& theImg);
     167             :   //casacore::Int storeTLAsImg(casacore::String fileName, casacore::TempLattice<casacore::Complex> &TL, casacore::ImageInterface<casacore::Float>& theImg);
     168             : 
     169             :   casacore::Bool mergeDataError(casacore::ImageInterface<casacore::Float> &data, casacore::ImageInterface<casacore::Float> &error, const casacore::String &outImg);
     170             : 
     171             :   casacore::Bool createMask(casacore::LatticeExpr<casacore::Bool> &lemask, casacore::ImageInterface<casacore::Float> &outimage);
     172             : 
     173             :   casacore::Bool resizeWorkArrays(casacore::Int length);
     174             :   
     175             :   casacore::Int makeSpectralPSFs(SkyEquation& se, casacore::Bool writeToDisk);
     176             :    //casacore::Int addTo(casacore::Lattice<casacore::Float>& to, const casacore::Lattice<casacore::Float>& add, casacore::Float multiplier);
     177             :   casacore::Int writeResultsToDisk();
     178             :   casacore::Float computeFluxLimit(casacore::Float &fractionOfPsf);
     179             : 
     180             :   void blankOverlappingModels();
     181             :   void restoreOverlappingModels();
     182             : 
     183             :   void saveCurrentModels();
     184             :   
     185             :   casacore::Timer tmr1,tmr2;
     186             :   casacore::Int adbg;
     187             :   casacore::Int tdbg;
     188             :   casacore::Int ddbg;
     189             :   // Put in some progress metre here...
     190             :   
     191             :   
     192             : };
     193             : 
     194             : 
     195             : } //# NAMESPACE CASA - END
     196             : 
     197             : #endif
     198             : 
     199             :