LCOV - code coverage report
Current view: top level - synthesis/ImagerObjects - CubeMinorCycleAlgorithm.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 177 0.0 %
Date: 2024-10-29 13:38:20 Functions: 0 11 0.0 %

          Line data    Source code
       1             : //# CubeMinorCycleAlgorithm.cc: implementation of class to deconvolve cube in parallel/serial 
       2             : //# Copyright (C) 2019
       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 General Public License as published by
       7             : //# the Free Software Foundation; either version 3 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 General Public
      13             : //# License for more details.
      14             : //#
      15             : //# https://www.gnu.org/licenses/
      16             : //#
      17             : //# Queries concerning CASA should be submitted at
      18             : //#        https://help.nrao.edu
      19             : //#
      20             : //#        Postal address: CASA Project Manager 
      21             : //#                        National Radio Astronomy Observatory
      22             : //#                        520 Edgemont Road
      23             : //#                        Charlottesville, VA 22903-2475 USA
      24             : //#
      25             : //#
      26             : //# $Id$
      27             : #include <casacore/lattices/Lattices/LatticeLocker.h>
      28             : #include <synthesis/ImagerObjects/CubeMinorCycleAlgorithm.h>
      29             : #include <synthesis/ImagerObjects/SynthesisDeconvolver.h>
      30             : #include <casacore/casa/Containers/Record.h>
      31             : #include <synthesis/ImagerObjects/SimpleSIImageStore.h>
      32             : #include <imageanalysis/Utilities/SpectralImageUtil.h>
      33             : 
      34             : using namespace casacore;
      35             : namespace casa { //# NAMESPACE CASA - BEGIN
      36             : extern Applicator applicator;
      37             : 
      38           0 :   CubeMinorCycleAlgorithm::CubeMinorCycleAlgorithm() : myName_p("CubeMinorCycleAlgorithm"), autoMaskOn_p(False),chanFlag_p(0), status_p(False){
      39             :         
      40           0 : }
      41           0 : CubeMinorCycleAlgorithm::~CubeMinorCycleAlgorithm() {
      42             :         
      43           0 : }
      44             :         
      45           0 : void CubeMinorCycleAlgorithm::get() {
      46             :   ///New instructions reset previous state
      47           0 :   reset();
      48             :         //cerr << "in get for child process " << applicator.isWorker() << endl;
      49           0 :         Record decParsRec;
      50             :      
      51             :         // get deconv params record #1
      52           0 :         applicator.get(decParsRec);
      53             :         // get iter control rec #2
      54           0 :         applicator.get(iterBotRec_p);
      55             :         // channel range to deconvolve #3
      56           0 :         applicator.get(chanRange_p);
      57             :         //get psf image name #4
      58           0 :         applicator.get(psfName_p);
      59             :         //get residual name #5
      60           0 :         applicator.get(residualName_p);
      61             :         //get model name #6
      62           0 :         applicator.get(modelName_p);
      63             :         //get mask name #7
      64           0 :         applicator.get(maskName_p);
      65             :         // get pb name #8
      66           0 :         applicator.get(pbName_p);
      67             :         //get beamsetrec #9
      68             :         //applicator.get(beamsetRec_p);
      69             :         //get psfsidelobelev #9
      70           0 :         applicator.get(psfSidelobeLevel_p);
      71             :         //get chanflag #10
      72           0 :         Record chanflagRec;
      73           0 :         applicator.get(chanflagRec);
      74           0 :         chanFlag_p.resize();
      75           0 :         chanflagRec.get("chanflag", chanFlag_p);
      76           0 :         statsRec_p=chanflagRec.asRecord("statsrec");
      77             :         //cerr <<"GET statsRec " << statsRec_p << endl;
      78           0 :         decPars_p.fromRecord(decParsRec);
      79             :         
      80             : 
      81             :         
      82           0 : }
      83           0 : void CubeMinorCycleAlgorithm::put() {
      84             :         
      85             :   ///# 1  chanrange processed 
      86           0 :   applicator.put(chanRange_p);
      87             :         //cerr << "in put " << status_p << endl;
      88             :   //#2 chanflag
      89           0 :   chanFlagRec_p.define("chanflag", chanFlag_p);
      90             :   //cerr << "PUT statsRec "<< statsRec_p << endl; 
      91           0 :   chanFlagRec_p.defineRecord("statsrec", statsRec_p);
      92           0 :   applicator.put(chanFlagRec_p);
      93             :   ///#3 return record of deconvolver
      94             :   // cerr << "nfield " << returnRec_p.nfields() << endl;
      95           0 :   applicator.put(returnRec_p);  
      96             :         
      97           0 : }
      98             :         
      99           0 : void CubeMinorCycleAlgorithm::task(){
     100           0 :         LogIO logger(LogOrigin("CubeMinorCycleAlgorithm", "task", WHERE));
     101             : 
     102           0 :         status_p = False;
     103             :         try{
     104             :          
     105           0 :           SynthesisDeconvolver subDeconv;
     106           0 :           Bool writeBackAutomask=True;
     107           0 :           subDeconv.setupDeconvolution(decPars_p);
     108           0 :           std::shared_ptr<SIImageStore> subimstor=subImageStore();
     109             :           //ImageBeamSet bs=ImageBeamSet::fromRecord(beamsetRec_p);
     110           0 :           ImageBeamSet bs=(subimstor->psf()->imageInfo()).getBeamSet();
     111           0 :           subimstor->setBeamSet(bs);
     112           0 :           subimstor->setPSFSidelobeLevel(psfSidelobeLevel_p);
     113           0 :           LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write, 30);
     114           0 :           Record prevresrec=subDeconv.initMinorCycle(subimstor);
     115             : 
     116             :           // unused
     117             :           // Float prevPeakRes=prevresrec.asFloat("peakresidual");
     118           0 :           Bool doDeconv=True;
     119           0 :           if(autoMaskOn_p){
     120           0 :             subDeconv.setChanFlag(chanFlag_p);
     121           0 :             subDeconv.setRobustStats(statsRec_p);
     122           0 :             Int automaskflag=iterBotRec_p.asInt("onlyautomask");
     123           0 :             LogIO os( LogOrigin("SynthesisDeconvolver",automaskflag==1 ? "excuteAutomask" : "executeMinorCycle",WHERE) );
     124           0 :             os << "Processing channels in range " << chanRange_p << LogIO::POST;
     125           0 :             if(automaskflag==1){
     126           0 :               doDeconv=False;
     127           0 :               if(iterBotRec_p.isDefined("cycleniter"))
     128           0 :                  subDeconv.setMinorCycleControl(iterBotRec_p);
     129             :             }
     130             :             //cerr << "ITERDONE " << iterBotRec_p.asInt("iterdone")<< " itermask flag " << automaskflag << endl;
     131           0 :             subDeconv.setIterDone(iterBotRec_p.asInt("iterdone"));
     132           0 :             if(automaskflag !=0){
     133             :               //this is already sent in as part of subimstor
     134             :               //subDeconv.setPosMask(subimstor->tempworkimage());
     135             :              
     136           0 :               subDeconv.setAutoMask();
     137             :               /*
     138             :               Record resRec=subDeconv.initMinorCycle(subimstor);
     139             :               //cerr << "POST resrec " << resRec << endl;
     140             :               //cerr << "peakResidual " << resRec.asFloat("peakresidual") << " cyclethreshold " << iterBotRec_p.asFloat("cyclethreshold") << " as double " << iterBotRec_p.asDouble("cyclethreshold") << endl;
     141             :               if(resRec.isDefined("peakresidual") && iterBotRec_p.isDefined("cyclethreshold")){
     142             :                 
     143             :                 Float peakresidual=resRec.asFloat("peakresidual");
     144             :                 Float peakresidualnomask=resRec.asFloat("peakresidualnomask");
     145             :                 Float cyclethreshold=iterBotRec_p.asFloat("cyclethreshold");
     146             :                 if(iterBotRec_p.isDefined("psffraction")){
     147             :                   cyclethreshold=iterBotRec_p.asFloat("psffraction")*peakresidual;    
     148             :                   cyclethreshold= max(cyclethreshold, iterBotRec_p.asFloat("threshold"));
     149             :                   if(automaskflag==-1){
     150             : 
     151             :                     //              cerr << "old cyclethreshold " <<iterBotRec_p.asFloat("cyclethreshold") << " new " << cyclethreshold << endl;
     152             :                     iterBotRec_p.removeField("cyclethreshold");
     153             :                     iterBotRec_p.define("cyclethreshold", cyclethreshold);
     154             :                   }
     155             :                 }
     156             :                 
     157             :                 
     158             :                   //if(peakresidual < iterBotRec_p.asFloat("cyclethreshold"))
     159             :                 if(peakresidual < cyclethreshold)
     160             :                   writeBackAutomask=False;
     161             :                 
     162             :                 writeBackAutomask=False;
     163             :                 //Its better to always write the automask 
     164             :                 //cerr << "chanRange " << chanRange_p << endl;
     165             :                 //cerr << "WRITEBACK " << writeBackAutomask<< " peakres " <<  peakresidual << " prev " <<  prevPeakRes << " nomask " << peakresidualnomask << endl;
     166             :                 
     167             :               }
     168             :               */
     169             : 
     170           0 :               writeBackAutomask=True;
     171             :               //Its better to always write the automask 
     172             :              
     173             :             }
     174             :             else{
     175           0 :               writeBackAutomask=False;
     176             :             }
     177             :             
     178           0 :           } else {
     179           0 :             LogIO os( LogOrigin("CubeMinorCycleAlgorithm","task",WHERE) );
     180           0 :             os << "Processing channels in range " << chanRange_p << LogIO::POST;
     181           0 :           } // end if(autoMaskOn_p)
     182             :           //subDeconv.setupMask();
     183           0 :           if(doDeconv)
     184           0 :             returnRec_p=subDeconv.executeCoreMinorCycle(iterBotRec_p);
     185             :           else
     186           0 :             returnRec_p.define("doneautomask", True);
     187           0 :           chanFlag_p.resize();
     188           0 :           chanFlag_p=subDeconv.getChanFlag();
     189           0 :           statsRec_p=Record();
     190           0 :           statsRec_p=subDeconv.getRobustStats();
     191           0 :           if(doDeconv){
     192           0 :             writeBackToFullImage(modelName_p, chanRange_p[0], chanRange_p[1], (subimstor->model()));
     193           0 :             writeBackToFullImage(residualName_p, chanRange_p[0], chanRange_p[1], (subimstor->residual()));
     194             : 
     195             :           }
     196           0 :           if(autoMaskOn_p && writeBackAutomask){
     197           0 :             writeBackToFullImage(posMaskName_p, chanRange_p[0], chanRange_p[1], (subimstor->tempworkimage()));
     198           0 :             writeBackToFullImage(maskName_p, chanRange_p[0], chanRange_p[1], (subimstor->mask()));
     199             :           }
     200           0 :         }
     201           0 :         catch (AipsError& x) {
     202           0 :           logger << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
     203           0 :           returnRec_p=Record();
     204           0 :         }
     205           0 :         catch(std::exception& exc) {
     206           0 :           logger << LogIO::SEVERE << "Exception (std): " << exc.what() << LogIO::POST;
     207           0 :           returnRec_p=Record();
     208           0 :         }
     209           0 :         catch(...){
     210           0 :           logger << LogIO::SEVERE << "Unknown exception" << LogIO::POST;
     211           0 :           returnRec_p=Record();
     212           0 :         }
     213             :         
     214           0 :         status_p = True;
     215           0 : }
     216           0 : String&     CubeMinorCycleAlgorithm::name(){
     217           0 :         return myName_p;
     218             : }
     219             : 
     220           0 : std::shared_ptr<SIImageStore> CubeMinorCycleAlgorithm::subImageStore(){
     221           0 :   std::shared_ptr<ImageInterface<Float> >subpsf=nullptr;
     222           0 :   std::shared_ptr<ImageInterface<Float> >subresid=nullptr;
     223           0 :   std::shared_ptr<ImageInterface<Float> >submodel=nullptr;
     224           0 :   std::shared_ptr<ImageInterface<Float> > submask=nullptr;
     225           0 :   std::shared_ptr<ImageInterface<Float> > subpb=nullptr;
     226           0 :   std::shared_ptr<ImageInterface<Float> > subposmask=nullptr;
     227           0 :         Int chanBeg=0;
     228           0 :         Int chanEnd=0;
     229           0 :         chanBeg=chanRange_p[0];
     230           0 :         chanEnd=chanRange_p[1];
     231           0 :         casacore::String imageName = decPars_p.imageName;
     232             :         //cerr << "chanBeg " << chanBeg << " chanEnd " << chanEnd << " imId " << imId << endl;
     233             :         
     234             :         
     235           0 :         makeTempImage(subpsf, psfName_p, chanBeg, chanEnd);
     236             :         ///Have to lock residual now that deconvolver task now will want to writeback residual
     237             :         ///tclean does not need writeback for residual  so it is unnecessary step for now.
     238           0 :         makeTempImage(subresid, residualName_p, chanBeg, chanEnd, True);
     239           0 :         makeTempImage(submodel, modelName_p, chanBeg, chanEnd, True);
     240           0 :         makeTempImage(submask, maskName_p, chanBeg, chanEnd, True);
     241             :         //       PagedImage<Float> model(modelName_p, TableLock::UserLocking);
     242             :         //   submodel.reset(SpectralImageUtil::getChannel(model, chanBeg, chanEnd, true));
     243           0 :         if(!pbName_p.empty()){
     244           0 :            makeTempImage(subpb, pbName_p, chanBeg, chanEnd);
     245             :         }
     246           0 :         if(iterBotRec_p.isDefined("posmaskname") ){
     247           0 :             iterBotRec_p.get("posmaskname", posMaskName_p);
     248           0 :             if(Table::isReadable(posMaskName_p)){
     249           0 :                 makeTempImage(subposmask, posMaskName_p, chanBeg, chanEnd, True);
     250           0 :                 if(subposmask)
     251           0 :                   autoMaskOn_p=True;
     252             :               }
     253             :           }
     254             : 
     255           0 :             std::shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imageName, submodel, subresid, subpsf, nullptr, nullptr, submask, nullptr, nullptr, subpb, nullptr, subposmask));
     256             :         
     257             :         //cerr << "subimagestor TYPE" << subimstor->getType() << endl;
     258           0 :         return subimstor;
     259           0 : }
     260           0 :   void CubeMinorCycleAlgorithm::makeTempImage(std::shared_ptr<ImageInterface<Float> >& outptr,  const String& imagename, const Int chanBeg, const Int chanEnd, const Bool writelock){
     261             :     //For testing that locks are placed in the right places use userlocking
     262           0 :     PagedImage<Float> im(imagename, writelock ? TableLock::UserLocking : TableLock::UserNoReadLocking);
     263             :     //PagedImage<Float> im(imagename, writelock ? TableLock::AutoLocking : TableLock::AutoNoReadLocking);
     264             :     
     265           0 :     SubImage<Float> *tmpptr=nullptr;
     266             :     //LatticeLocker lockread (im, FileLocker::Read);
     267             :     //if(writelock)
     268             :     //  im.lock(FileLocker::Write, 1000);
     269             :     ////TESTOO
     270             :     //outptr.reset(SpectralImageUtil::getChannel(im, chanBeg, chanEnd, writelock));
     271             :     
     272             : 
     273             :     ///END of TESTOO
     274           0 :     if(writelock){
     275           0 :       im.lock(FileLocker::Write, 1000);
     276           0 :       tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
     277             :     }
     278             :     else
     279           0 :       tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
     280           0 :     if(tmpptr){
     281           0 :       IPosition tileshape=tmpptr->shape();
     282           0 :       tileshape[2]=1; tileshape[3]=1;
     283           0 :       TiledShape tshape(tmpptr->shape(),tileshape);
     284             :       ///TESTOO
     285             :       /*if(imagename.contains(".residual")){
     286             :         String newresid=File::newUniqueName("./", "newResid").absoluteName();
     287             :         outptr.reset(new PagedImage<Float>(tshape, tmpptr->coordinates(), newresid));
     288             :         }
     289             :       ////TESTOO
     290             :       else
     291             :       */
     292           0 :       outptr.reset(new TempImage<Float>(tshape, tmpptr->coordinates()));
     293             :         
     294             :       
     295           0 :       if(writelock){
     296           0 :         LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
     297           0 :         outptr->copyData(*tmpptr);
     298             :         //cerr << "IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr  " << tmpptr->isMasked() << endl;
     299           0 :         if(tmpptr->isMasked()){
     300           0 :           outptr->makeMask ("mask0", true, true, false, true);
     301           0 :           outptr->pixelMask().put(tmpptr->getMask());
     302             :           //      cerr << "tempimage SUM of bit mask" << ntrue(tmpptr->pixelMask().get()) << " out   " << ntrue(outptr->pixelMask().get()) << endl;
     303             :         }
     304           0 :       }
     305             :       else{
     306           0 :         LatticeLocker lock1 (*(tmpptr), FileLocker::Read);
     307           0 :         outptr->copyData(*tmpptr);
     308             :         //cerr << "false IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr  " << tmpptr->isMasked() << endl;
     309             :         
     310           0 :         if(tmpptr->isMasked()){
     311           0 :           outptr->makeMask ("mask0", true, true, false, true);
     312           0 :           outptr->pixelMask().put(tmpptr->getMask());
     313             :           //cerr << "tempimage SUM of bit mask" << ntrue(tmpptr->pixelMask().get()) << " out   " << ntrue(outptr->pixelMask().get()) << endl;
     314             :         }
     315           0 :       }
     316           0 :       ImageInfo iinfo=tmpptr->imageInfo();
     317           0 :       outptr->setImageInfo(iinfo);
     318           0 :       delete tmpptr;
     319           0 :     }
     320             :     
     321           0 :     im.unlock();
     322           0 :   }
     323           0 :  void CubeMinorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
     324           0 :    PagedImage<Float> im(imagename, TableLock::UserLocking);
     325             :    //PagedImage<Float> im(imagename, TableLock::AutoLocking);
     326           0 :     SubImage<Float> *tmpptr=nullptr; 
     327           0 :     tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, true);
     328             :     {
     329             :     
     330           0 :       LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
     331           0 :       tmpptr->copyData(*(subimptr));
     332           0 :     }        
     333           0 :     im.unlock();
     334           0 :     delete tmpptr;
     335             :                  
     336           0 :   }
     337           0 : void CubeMinorCycleAlgorithm::reset(){
     338             :                 
     339           0 :   iterBotRec_p=Record();
     340           0 :   modelName_p="";
     341           0 :   residualName_p="";
     342           0 :   psfName_p="";
     343           0 :   maskName_p="";
     344           0 :   pbName_p="";
     345           0 :   posMaskName_p="";
     346           0 :   chanRange_p.resize();
     347           0 :   returnRec_p=Record();
     348           0 :   beamsetRec_p=Record();
     349             :   //psfSidelobeLevel_p;
     350           0 :   autoMaskOn_p=False;
     351           0 :   chanFlag_p.resize();
     352           0 :   statsRec_p=Record();
     353           0 :   chanFlagRec_p=Record();
     354           0 :   status_p=False;
     355             :                 
     356             :         
     357             :         
     358           0 : }       
     359             :         
     360             :         
     361             : } //# NAMESPACE CASA - END

Generated by: LCOV version 1.16