Line data    Source code

       1             : //# ClarkCleanImageSkyModel.cc: Implementation of ClarkCleanImageSkyModel class
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
       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: casa-feedback@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             : #include <casacore/casa/Arrays/ArrayMath.h>
      29             : #include <synthesis/MeasurementComponents/ClarkCleanImageSkyModel.h>
      30             : #include <casacore/images/Images/PagedImage.h>
      31             : #include <casacore/casa/OS/File.h>
      32             : #include <casacore/lattices/Lattices/LatticeStepper.h>
      33             : #include <casacore/lattices/Lattices/LatticeIterator.h>
      34             : #include <casacore/lattices/LEL/LatticeExpr.h>
      35             : #include <casacore/lattices/LEL/LatticeExprNode.h>
      36             : #include <casacore/lattices/LRegions/LCBox.h>
      37             : #include <casacore/lattices/Lattices/SubLattice.h>
      38             : #include <synthesis/MeasurementEquations/SkyEquation.h>
      39             : #include <casacore/casa/Exceptions/Error.h>
      40             : #include <casacore/casa/BasicSL/String.h>
      41             : #include <casacore/casa/Utilities/Assert.h>
      42             : 
      43             : #include <sstream>
      44             : 
      45             : #include <casacore/casa/Logging/LogMessage.h>
      46             : #include <casacore/casa/Logging/LogIO.h>
      47             : #include <casacore/casa/Logging/LogSink.h>
      48             : 
      49             : #include <synthesis/MeasurementEquations/LatConvEquation.h>
      50             : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
      51             : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
      52             : 
      53             : 
      54             : using namespace casacore;
      55             : namespace casa { //# NAMESPACE CASA - BEGIN
      56             : 
      57           0 :   ClarkCleanImageSkyModel::ClarkCleanImageSkyModel() : CleanImageSkyModel(), itsProgress(0) {
      58             : 
      59             : 
      60             : 
      61           0 : }
      62             : 
      63             : 
      64           0 : ClarkCleanImageSkyModel::~ClarkCleanImageSkyModel()
      65             : { 
      66           0 :   if (itsProgress) delete itsProgress; 
      67           0 : };
      68             : 
      69             : // Clean solver
      70           0 : Bool ClarkCleanImageSkyModel::solve(SkyEquation& se) {
      71             : 
      72           0 :   LogIO os(LogOrigin("ClarkCleanImageSkyModel","solve",WHERE));
      73             :   
      74           0 :   Bool converged=false; 
      75           0 :   if(numberOfModels()>1) {
      76           0 :     os << "Cannot process more than one field" << LogIO::EXCEPTION;
      77             :   }
      78             : 
      79             :   // Make the residual image
      80           0 :   if(modified_p)
      81           0 :     makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
      82             :   
      83             :   //Make the PSF
      84           0 :   if(!donePSF_p)
      85           0 :     makeApproxPSFs(se);
      86             : 
      87           0 :   if(numberIterations() < 1){
      88             : 
      89           0 :     return true;
      90             :   }
      91             : 
      92           0 :   Float maxRes=0.0;
      93           0 :   Int iterused=0;
      94             :   
      95           0 :   if(hasMask(0))
      96           0 :     converged = clean(image(0), residual(0), PSF(0), mask(0), maxRes, iterused, gain(), numberIterations(),  
      97           0 :           threshold(), cycleFactor_p, true, doPolJoint_p);
      98             :   else{
      99           0 :     ImageInterface<Float> *tmpMask=0;
     100           0 :     converged = clean(image(0), residual(0), PSF(0), *tmpMask, maxRes, iterused, gain(), numberIterations(),  
     101           0 :           threshold(), cycleFactor_p, false, doPolJoint_p);
     102             :   }
     103             :   
     104           0 :   setThreshold(maxRes);
     105           0 :   setNumberIterations(iterused);
     106           0 :   modified_p=true;
     107           0 :   return(converged);
     108           0 : };
     109             : 
     110             : 
     111           0 :   Bool ClarkCleanImageSkyModel::clean(ImageInterface<Float>& image, ImageInterface<Float> & residual, 
     112             :                                       ImageInterface<Float>& psf, ImageInterface<Float>& mask, Float& maxresidual, Int& iterused, Float gain, Int numIter,  
     113             :                                       Float thresh, Float cycleFactor, Bool useMask, Bool doPolJoint){
     114             : 
     115           0 :     Bool converged=false;
     116           0 :     LogIO os(LogOrigin("ClarkCleanImageSkyModel","clean",WHERE)); 
     117           0 :     Int nx=image.shape()(0);
     118           0 :     Int ny=image.shape()(1);
     119           0 :     Int npol=image.shape()(2);
     120           0 :     Int nchan=image.shape()(3);
     121           0 :     Int polloop=1;
     122           0 :     Vector<String> stokesID(npol, "");
     123           0 :     if (!doPolJoint) {
     124           0 :       polloop=npol;
     125           0 :       CoordinateSystem cs=image.coordinates();
     126           0 :       Int stokesindex=cs.findCoordinate(Coordinate::STOKES);
     127           0 :       StokesCoordinate stcoord=cs.stokesCoordinate(stokesindex);
     128           0 :       for (Int jj=0; jj < npol ; ++jj){
     129           0 :         stokesID(jj)=Stokes::name(Stokes::type(stcoord.stokes()(jj)));      
     130             :       }
     131           0 :     }
     132             :   
     133             :   
     134           0 :     Int newNx=nx;
     135           0 :     Int newNy=ny;
     136             : 
     137             :     Int xbeg, xend, ybeg, yend;
     138             :     //default clean box
     139           0 :     xbeg=nx/4; 
     140           0 :     xend=3*nx/4-1;
     141           0 :     ybeg=ny/4; 
     142           0 :     yend=3*ny/4-1;
     143             : 
     144           0 :     Bool isCubeMask=false; 
     145             :     //AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
     146             :   
     147           0 :     Lattice<Float>* mask_sl = 0;
     148           0 :     RO_LatticeIterator<Float>* maskli = 0;
     149             :   
     150           0 :     if(useMask) {
     151             :       // AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
     152             :       // AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
     153           0 :       if((mask.shape()(0)!=nx) || (mask.shape()(1)!=ny)){
     154           0 :         throw(AipsError("Mask image shape is different from dirty image"));
     155             :       }
     156           0 :       if(nchan >1){
     157           0 :         if(mask.shape()(3)==nchan){
     158           0 :           isCubeMask=true;
     159           0 :           os << "Using multichannel mask" << LogIO::POST;
     160             :         }
     161             :         else{
     162             :           os << "Image cube and mask donot match in number of channels" 
     163           0 :              << LogIO::WARN;
     164             :           os << "Will use first plane of the mask for all channels" 
     165           0 :              << LogIO::WARN;
     166             :         }
     167             :       }
     168           0 :       LatticeStepper mls(mask.shape(),
     169           0 :                          IPosition(4, nx, ny, 1, 1),
     170           0 :                          IPosition(4, 0, 1, 3, 2));
     171           0 :       maskli= new RO_LatticeIterator<Float>(mask, mls);
     172           0 :       maskli->reset();
     173           0 :       mask_sl=makeMaskSubLat(nx, ny, newNx, newNy, *maskli, xbeg, xend, 
     174             :                              ybeg, yend);
     175           0 :     }
     176             : 
     177             : 
     178           0 :     maxresidual=0.0;
     179             : 
     180           0 :     for (Int ichan=0; ichan < nchan; ichan++) {
     181           0 :       if(useMask && isCubeMask && ichan >0) {
     182           0 :         (*maskli)++;
     183           0 :         if(mask_sl) delete mask_sl;
     184           0 :         mask_sl=0;
     185           0 :         mask_sl=makeMaskSubLat(nx, ny, newNx, newNy, *maskli, xbeg, 
     186             :                                xend, ybeg, yend);       
     187             :       }
     188             : 
     189           0 :       if(nchan>1) {
     190           0 :         os<<"Processing channel "<<ichan+1<<" of "<<nchan<<LogIO::POST;
     191             :       }
     192           0 :       for (Int ipol=0; ipol < polloop; ++ipol){
     193           0 :         Int polbeg=0;
     194           0 :         Int polend=npol-1;
     195           0 :         if(!doPolJoint){
     196           0 :           polbeg=ipol;
     197           0 :           polend=ipol;
     198           0 :           os << "Doing stokes "<< stokesID(ipol) << " image" <<LogIO::POST;
     199             :         }
     200           0 :         LCBox imagebox(IPosition(4, xbeg, ybeg, polbeg, ichan), 
     201           0 :                        IPosition(4, xend, yend, polend, ichan),
     202           0 :                        image.shape());
     203           0 :         LCBox psfbox(IPosition(4, 0, 0, 0, ichan), 
     204           0 :                      IPosition(4, nx-1, ny-1, 0, ichan),
     205           0 :                      psf.shape());
     206             : 
     207           0 :         SubLattice<Float> psf_sl(psf, psfbox, false);
     208           0 :         SubLattice<Float>  residual_sl(residual, imagebox, true);
     209           0 :         SubLattice<Float>  model_sl=SubLattice<Float>   (image, imagebox, true);
     210             :     
     211             :     
     212           0 :         ArrayLattice<Float> psftmp;
     213             :         
     214           0 :         if(nx != newNx){
     215             :           
     216           0 :           psftmp=ArrayLattice<Float> (IPosition(4, newNx, newNy, 1, 1));
     217           0 :           psftmp.set(0.0);
     218           0 :           Array<Float> tmparr=psf_sl.get();
     219           0 :           psftmp.putSlice(tmparr, IPosition(4, (newNx-nx)/2, (newNy-ny)/2, 0,0));
     220           0 :           psf_sl=SubLattice<Float>(psftmp, true);
     221             :  
     222           0 :         }
     223             :       
     224           0 :         TempLattice<Float> dirty_sl( residual_sl.shape());
     225           0 :         dirty_sl.copyData(residual_sl);
     226           0 :         TempLattice<Float> localmodel(model_sl.shape());
     227           0 :         localmodel.set(0.0);
     228             :       
     229             :         Float psfmax;
     230             :         {
     231           0 :           LatticeExprNode node = max(psf_sl);
     232           0 :           psfmax = node.getFloat();
     233           0 :         }
     234             :       
     235           0 :         if((psfmax==0.0) ||(useMask && (mask_sl == 0)) ) {
     236           0 :           maxresidual=0.0;
     237           0 :           iterused=0;
     238             :           os << LogIO::NORMAL // Loglevel INFO
     239           0 :              << "No data or blank mask for this channel: skipping" << LogIO::POST;
     240             :         } else {
     241           0 :           LatConvEquation eqn(psf_sl, residual_sl);
     242           0 :           ClarkCleanLatModel cleaner( localmodel );
     243           0 :           cleaner.setResidual(dirty_sl);
     244           0 :           if (mask_sl != 0 ) cleaner.setMask( *mask_sl );
     245             :         
     246             :           /*
     247             :           ClarkCleanProgress *cpp  = 0;
     248             :           if (displayProgress_p) {
     249             :             cpp = new ClarkCleanProgress( pgplotter_p );
     250             :             cleaner.setProgress(*cpp);
     251             :           }
     252             :           */
     253             :         
     254           0 :           cleaner.setGain(gain);
     255           0 :           cleaner.setNumberIterations(numIter);
     256           0 :           cleaner.setThreshold(thresh); 
     257           0 :           cleaner.setPsfPatchSize(IPosition(2,51,51)); 
     258           0 :           cleaner.setHistLength(1024);
     259           0 :           cleaner.setMaxNumPix(32*1024);
     260           0 :           cleaner.setCycleFactor(cycleFactor);
     261             :           // clean if there is no mask or if it has mask AND mask is not empty 
     262           0 :           cleaner.solve(eqn);
     263           0 :           cleaner.setChoose(false);
     264           0 :           if (cleaner.numberIterations()>0) {
     265             :             os << LogIO::NORMAL // Loglevel INFO
     266             :                << "Clean used " << cleaner.numberIterations() << " iterations" 
     267           0 :                << " in this round to get to a max residual of " << cleaner.threshold()
     268           0 :                << LogIO::POST;
     269             :           } 
     270           0 :           LatticeExpr<Float> expr= model_sl + localmodel; 
     271           0 :           model_sl.copyData(expr);
     272           0 :           maxresidual=max(maxresidual, cleaner.getMaxResidual());
     273           0 :           iterused=cleaner.numberIterations();
     274           0 :           converged =  (cleaner.getMaxResidual() < thresh) 
     275           0 :             || (cleaner.numberIterations()==0);
     276             :           //      if (cpp != 0 ) delete cpp; cpp=0;
     277             :           //      if (pgp != 0 ) delete pgp; pgp=0;
     278           0 :         }
     279           0 :       }
     280             :     }
     281           0 :     if (mask_sl != 0)  {
     282           0 :       delete mask_sl;
     283           0 :       mask_sl=0;
     284             :     }
     285           0 :     if (maskli !=0) {
     286           0 :       delete maskli;
     287           0 :       maskli=0;
     288             :     }
     289             :     os << LogIO::NORMAL // Loglevel INFO
     290           0 :        << LatticeExprNode(sum(image)).getFloat()
     291             :        << " Jy <- The sum of the clean components"
     292           0 :        << LogIO::POST;
     293           0 :     return converged;
     294             : 
     295             : 
     296           0 :   }
     297             : 
     298           0 : Lattice<Float>* ClarkCleanImageSkyModel::makeMaskSubLat(const Int& nx, 
     299             :                                                             const Int& ny, 
     300             :                                                         Int& newNx, Int& newNy,
     301             :                                                             RO_LatticeIterator<Float>& maskIter,
     302             :                                                            Int& xbeg,
     303             :                                                            Int& xend,
     304             :                                                            Int& ybeg,
     305             :                                                            Int& yend) {
     306             : 
     307           0 :   LogIO os(LogOrigin("ClarkCleanImageSkyModel","makeMaskSubLat",WHERE)); 
     308             : 
     309             : 
     310           0 :   newNx=nx;
     311           0 :   newNy=ny;
     312           0 :   SubLattice<Float>* mask_sl = 0;
     313           0 :   xbeg=nx/4;
     314           0 :   ybeg=ny/4;
     315             :   
     316           0 :   xend=xbeg+nx/2-1;
     317           0 :   yend=ybeg+ny/2-1;  
     318           0 :   Matrix<Float> mask= maskIter.matrixCursor();
     319           0 :   if(max(mask) < 0.000001) {
     320             :     //zero mask ....
     321           0 :     return mask_sl;
     322             :   }
     323             :   // Now read the mask and determine the bounding box
     324             : 
     325           0 :   xbeg=nx-1;
     326           0 :   ybeg=ny-1;
     327           0 :   xend=0;
     328           0 :   yend=0;
     329             : 
     330             :   
     331           0 :   for (Int iy=0;iy<ny;iy++) {
     332           0 :     for (Int ix=0;ix<nx;ix++) {
     333           0 :       if(mask(ix,iy)>0.000001) {
     334           0 :         xbeg=min(xbeg,ix);
     335           0 :         ybeg=min(ybeg,iy);
     336           0 :         xend=max(xend,ix);
     337           0 :         yend=max(yend,iy);
     338             : 
     339             :       }
     340             :     }
     341             :   }
     342             :   // Now have possible BLC. Make sure that we don't go over the
     343             :   // edge later
     344           0 :   Bool larger_quarter=false;
     345           0 :   if((xend - xbeg)>nx/2) {
     346             :     // xbeg=nx/4-1; //if larger than quarter take inner of mask
     347             :     //os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis"  << LogIO::POST;
     348           0 :     larger_quarter=true;
     349             :   } 
     350           0 :   if((yend - ybeg)>ny/2) { 
     351             :     //ybeg=ny/4-1;
     352             :     //os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
     353           0 :     larger_quarter=true;
     354             :   }  
     355             : 
     356             :   // make sure xend, yend does not go out of bounds
     357           0 :   if(larger_quarter){
     358           0 :     xend=min(xend, nx-1);
     359           0 :     yend=min(yend, ny-1);
     360             : 
     361             :   }
     362             :   else{
     363           0 :     xend=min(xend,xbeg+nx/2-1);
     364           0 :     yend=min(yend,ybeg+ny/2-1); 
     365             :   }
     366             : 
     367             :   
     368           0 :   if ((xend > xbeg) && (yend > ybeg) ) {
     369             : 
     370           0 :     IPosition latshape(4, mask.shape()(0), mask.shape()(1), 1,1);
     371           0 :     ArrayLattice<Float> arrayLat(latshape);
     372           0 :     LCBox maskbox (IPosition(4, xbeg, ybeg, 0, 0), 
     373           0 :                    IPosition(4, xend, yend, 0, 0), 
     374           0 :                    latshape);
     375             :     
     376           0 :     arrayLat.putSlice(mask, IPosition(4,0,0,0,0));
     377           0 :     mask_sl=new SubLattice<Float> (arrayLat, maskbox, false);
     378             :     
     379             : 
     380             : 
     381           0 :   }
     382             :   
     383           0 :   if(larger_quarter){
     384           0 :     newNx=2*nx;
     385           0 :     newNy=2*ny;
     386             :     /*xbeg=xbeg+(newNx-nx)/2;
     387             :     ybeg=ybeg+(newNy-ny)/2;
     388             :     xend=xend+(newNx-nx)/2;
     389             :     yend=yend+(newNy-ny)/2;
     390             :     */
     391             :   }
     392           0 :   return mask_sl;
     393           0 : }
     394             : 
     395             : 
     396             : 
     397             : 
     398             : 
     399             : } //# NAMESPACE CASA - END
     400             :