Line data Source code
1 : //# SynthesisDeconvolver.cc: Implementation of Imager.h
2 : //# Copyright (C) 1997-2020
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 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/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 :
48 : #include <casacore/images/Images/TempImage.h>
49 : #include <casacore/images/Images/SubImage.h>
50 : #include <casacore/images/Regions/ImageRegion.h>
51 : #include <casacore/lattices/Lattices/LatticeLocker.h>
52 : #include <synthesis/ImagerObjects/CubeMinorCycleAlgorithm.h>
53 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
54 : #include <synthesis/ImagerObjects/SynthesisDeconvolver.h>
55 : #include <synthesis/ImagerObjects/SIMinorCycleController.h>
56 :
57 : #include <sys/types.h>
58 : #include <unistd.h>
59 : using namespace std;
60 :
61 : using namespace casacore;
62 : extern casa::Applicator casa::applicator;
63 : namespace casa { //# NAMESPACE CASA - BEGIN
64 :
65 17 : SynthesisDeconvolver::SynthesisDeconvolver() :
66 17 : itsDeconvolver(nullptr),
67 17 : itsMaskHandler(nullptr ),
68 17 : itsImages(nullptr),
69 17 : itsImageName(""),
70 : //itsPartImageNames(Vector<String>(0)),
71 17 : itsBeam(0.0),
72 17 : itsDeconvolverId(0),
73 17 : itsScales(Vector<Float>()),
74 17 : itsMaskType(""),
75 17 : itsPBMask(0.0),
76 : //itsMaskString(String("")),
77 17 : itsIterDone(0.0),
78 17 : itsChanFlag(Vector<Bool>(0)),
79 17 : itsRobustStats(Record()),
80 17 : initializeChanMaskFlag(false),
81 17 : itsPosMask(nullptr),
82 17 : itsIsMaskLoaded(false),
83 17 : itsMaskSum(-1e+9),
84 17 : itsPreviousFutureRes(0.0),
85 85 : itsPreviousIterBotRec_p(Record())
86 : {
87 17 : }
88 :
89 17 : SynthesisDeconvolver::~SynthesisDeconvolver()
90 : {
91 34 : LogIO os( LogOrigin("SynthesisDeconvolver","descructor",WHERE) );
92 17 : os << LogIO::DEBUG1 << "SynthesisDeconvolver destroyed" << LogIO::POST;
93 17 : SynthesisUtilMethods::getResource("End SynthesisDeconvolver");
94 :
95 17 : }
96 :
97 17 : void SynthesisDeconvolver::setupDeconvolution(const SynthesisParamsDeconv& decpars)
98 : {
99 34 : LogIO os( LogOrigin("SynthesisDeconvolver","setupDeconvolution",WHERE) );
100 :
101 : //Copy this decpars into a private variable that can be used elsewhere
102 : //there is no proper copy operator (as public casa::Arrays members = operator fails)
103 17 : itsDecPars.fromRecord(decpars.toRecord());
104 17 : itsImageName = decpars.imageName;
105 17 : itsStartingModelNames = decpars.startModel;
106 17 : itsDeconvolverId = decpars.deconvolverId;
107 :
108 17 : os << "Set Deconvolution Options for [" << itsImageName << "] : " << decpars.algorithm ;
109 :
110 17 : if( itsStartingModelNames.nelements()>0 && itsStartingModelNames[0].length() > 0 )
111 0 : os << " , starting from model : " << itsStartingModelNames;
112 17 : os << LogIO::POST;
113 :
114 : try
115 : {
116 17 : if(decpars.algorithm==String("hogbom"))
117 : {
118 0 : itsDeconvolver.reset(new SDAlgorithmHogbomClean());
119 : }
120 17 : else if(decpars.algorithm==String("mtmfs"))
121 : {
122 0 : itsDeconvolver.reset(new SDAlgorithmMSMFS( decpars.nTaylorTerms, decpars.scales, decpars.scalebias ));
123 : }
124 17 : else if(decpars.algorithm==String("clark_exp"))
125 : {
126 0 : itsDeconvolver.reset(new SDAlgorithmClarkClean("clark"));
127 : }
128 17 : else if(decpars.algorithm==String("clarkstokes_exp"))
129 : {
130 0 : itsDeconvolver.reset(new SDAlgorithmClarkClean("clarkstokes"));
131 : }
132 17 : else if(decpars.algorithm==String("clark"))
133 : {
134 17 : itsDeconvolver.reset(new SDAlgorithmClarkClean2("clark"));
135 : }
136 0 : else if(decpars.algorithm==String("clarkstokes"))
137 : {
138 0 : itsDeconvolver.reset(new SDAlgorithmClarkClean2("clarkstokes"));
139 : }
140 0 : else if(decpars.algorithm==String("multiscale"))
141 : {
142 0 : itsDeconvolver.reset(new SDAlgorithmMSClean( decpars.scales, decpars.scalebias ));
143 : }
144 0 : else if(decpars.algorithm==String("mem"))
145 : {
146 0 : itsDeconvolver.reset(new SDAlgorithmMEM( "entropy" ));
147 : }
148 0 : else if (decpars.algorithm==String("asp"))
149 : {
150 : /*bool isSingle = false;
151 : if (decpars.specmode == String("mfs"))
152 : isSingle = true;
153 : */
154 :
155 : // CAS-13901 : When `deconvolve()` is called on a cube,
156 : // `decpars.specmode` defaults to MFS here, since there is no way to
157 : // specify specmode within deconvolve. Therefore `isSingle` is
158 : // forced to False always. The side-effect of this is that
159 : // multi-plane MFS images will be slightly slowed down.
160 :
161 0 : bool isSingle = false;
162 0 : itsDeconvolver.reset(new SDAlgorithmAAspClean(decpars.fusedThreshold, isSingle, decpars.largestscale));
163 : }
164 : else
165 : {
166 0 : throw( AipsError("Un-known algorithm : "+decpars.algorithm) );
167 : }
168 :
169 : // Set restoring beam options
170 17 : itsDeconvolver->setRestoringBeam( decpars.restoringbeam, decpars.usebeam );
171 17 : itsUseBeam = decpars.usebeam;// store this info also here.
172 :
173 : // Set Masking options
174 : // itsDeconvolver->setMaskOptions( decpars.maskType );
175 17 : itsMaskHandler.reset(new SDMaskHandler());
176 : //itsMaskString = decpars.maskString;
177 17 : itsMaskType = decpars.maskType;
178 17 : if(itsMaskType=="auto-thresh")
179 : {
180 0 : itsAutoMaskAlgorithm="thresh";
181 : }
182 17 : else if(itsMaskType=="auto-thresh2")
183 : {
184 0 : itsAutoMaskAlgorithm="thresh2";
185 : }
186 17 : else if(itsMaskType=="auto-multithresh")
187 : {
188 0 : itsAutoMaskAlgorithm="multithresh";
189 : }
190 17 : else if(itsMaskType=="auto-onebox")
191 : {
192 0 : itsAutoMaskAlgorithm="onebox";
193 : }
194 : else {
195 17 : itsAutoMaskAlgorithm="";
196 : }
197 17 : itsPBMask = decpars.pbMask;
198 17 : itsMaskString = decpars.maskString;
199 51 : if(decpars.maskList.nelements()==0 ||
200 34 : (decpars.maskList.nelements()==1 && decpars.maskList[0]==""))
201 : {
202 17 : itsMaskList.resize(1);
203 17 : itsMaskList[0] = itsMaskString;
204 : }
205 : else
206 : {
207 0 : itsMaskList = decpars.maskList;
208 : }
209 17 : itsMaskThreshold = decpars.maskThreshold;
210 17 : itsFracOfPeak = decpars.fracOfPeak;
211 17 : itsMaskResolution = decpars.maskResolution;
212 17 : itsMaskResByBeam = decpars.maskResByBeam;
213 17 : itsNMask = decpars.nMask;
214 : //itsAutoAdjust = decpars.autoAdjust;
215 : //desable autoadjust
216 17 : itsAutoAdjust = false;
217 17 : itsSidelobeThreshold = decpars.sidelobeThreshold;
218 17 : itsNoiseThreshold = decpars.noiseThreshold;
219 17 : itsLowNoiseThreshold = decpars.lowNoiseThreshold;
220 17 : itsNegativeThreshold = decpars.negativeThreshold;
221 17 : itsSmoothFactor = decpars.smoothFactor;
222 17 : itsMinBeamFrac = decpars.minBeamFrac;
223 17 : itsCutThreshold = decpars.cutThreshold;
224 17 : itsGrowIterations = decpars.growIterations;
225 17 : itsDoGrowPrune = decpars.doGrowPrune;
226 17 : itsMinPercentChange = decpars.minPercentChange;
227 17 : itsVerbose = decpars.verbose;
228 17 : itsFastNoise = decpars.fastnoise;
229 17 : itsIsInteractive = decpars.interactive;
230 17 : itsNsigma = decpars.nsigma;
231 17 : itsNoRequireSumwt = true; //decpars.noRequireSumwt;
232 17 : itsFullSummary = decpars.fullsummary;
233 : }
234 0 : catch(AipsError &x)
235 : {
236 0 : throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
237 0 : }
238 :
239 17 : itsAddedModel=false;
240 17 : }
241 :
242 17 : Long SynthesisDeconvolver::estimateRAM(const vector<int>& imsize){
243 :
244 17 : Long mem=0;
245 : /* This does not work
246 : if( ! itsImages )
247 : {
248 : itsImages = makeImageStore( itsImageName );
249 : }
250 : */
251 17 : if(itsDeconvolver)
252 17 : mem=itsDeconvolver->estimateRAM(imsize);
253 17 : return mem;
254 : }
255 :
256 0 : Record SynthesisDeconvolver::initMinorCycle() {
257 : /////IMPORTANT initMinorCycle has to be called before setupMask...that order has to be kept !
258 :
259 0 : if(!itsImages)
260 0 : itsImages=makeImageStore(itsImageName, itsNoRequireSumwt);
261 : //For cubes as we are not doing a post major cycle residual automasking
262 : //Force recalculation of robust stats to update nsigmathreshold with
263 : //most recent residual
264 :
265 0 : if(itsAutoMaskAlgorithm=="multithresh" && itsImages->residual()->shape()[3] >1 && itsNsigma > 0.0){
266 0 : Record retval;
267 0 : Record backupRobustStats=itsRobustStats;
268 0 : itsRobustStats=Record();
269 0 : retval=initMinorCycle(itsImages);
270 0 : itsRobustStats=backupRobustStats;
271 0 : return retval;
272 0 : }
273 :
274 : /* else if (itsAutoMaskAlgorithm=="multithresh" && itsImages->residual()->shape()[3]){
275 : ///As the automask for cubes pre-CAS-9386...
276 : /// was tuned to look for threshold in future mask
277 : ///It is as good as somewhere in between no mask and mask
278 : // Record backupRobustStats=itsRobustStats;
279 : Record retval=initMinorCycle(itsImages);
280 : //cerr << "INITMINOR " << itsRobustStats << endl;
281 : //itsRobustStats=backupRobustStats;
282 : if(retval.isDefined("peakresidualnomask")){
283 : Float futureRes=Float(retval.asFloat("peakresidualnomask")-(retval.asFloat("peakresidualnomask")-retval.asFloat("peakresidual"))/1000.0);
284 : if(futureRes != itsPreviousFutureRes){
285 : //itsLoopController.setPeakResidual(retval.asFloat("peakresidualnomask"));
286 : retval.define("peakresidual", futureRes);
287 : itsPreviousFutureRes=futureRes;
288 : }
289 : }
290 : return retval;
291 : }
292 : */
293 0 : Record retval= initMinorCycle(itsImages);
294 : // cerr << "INITMINOR retval" << retval << endl;
295 :
296 0 : return retval;
297 0 : }
298 0 : Record SynthesisDeconvolver::initMinorCycle(std::shared_ptr<SIImageStore> imstor )
299 : {
300 0 : LogIO os( LogOrigin("SynthesisDeconvolver","initMinorCycle",WHERE) );
301 0 : Record returnRecord;
302 0 : Timer timer;
303 0 : Timer tim;
304 0 : tim.mark();
305 : try {
306 :
307 : //os << "---------------------------------------------------- Init (?) Minor Cycles ---------------------------------------------" << LogIO::POST;
308 :
309 : Int nSubChans, nSubPols;
310 0 : nSubChans = imstor->getShape()(3);
311 0 : nSubPols = imstor->getShape()(2);
312 0 : itsImages = imstor;
313 :
314 : // If a starting model exists, this will initialize the ImageStore with it. Will do this only once.
315 0 : setStartingModel();
316 :
317 : //itsIterDone is currently only used by automask code so move this to inside setAutomask
318 : //itsIterDone += itsLoopController.getIterDone();
319 :
320 : // setupMask();
321 : //
322 : Float masksum;
323 0 : if( ! itsImages->hasMask() ) // i.e. if there is no existing mask to re-use...
324 : {
325 0 : masksum = -1.0;
326 : }
327 : else
328 : {
329 0 : masksum = itsImages->getMaskSum();
330 0 : itsImages->mask()->unlock();
331 : }
332 0 : Bool validMask = ( masksum > 0 );
333 :
334 : // os << LogIO::NORMAL3 << "****INITMINOR Masksum stuff "<< tim.real() << LogIO::POST;
335 : // tim.mark();
336 :
337 : // Calculate Peak Residual and Max Psf Sidelobe, and fill into SubIterBot.
338 0 : Float peakresnomask = itsImages->getPeakResidual();
339 : // CAS-14201 : if the mask is all zero, then the peakresinmask is 0
340 0 : Float peakresinmask= validMask ? itsImages->getPeakResidualWithinMask() : 0;
341 : //os << LogIO::NORMAL3 << "****INITMINOR residual peak "<< tim.real() << LogIO::POST;
342 : //tim.mark();
343 0 : itsLoopController.setPeakResidual( validMask ? peakresinmask : peakresnomask );
344 : //os << LogIO::NORMAL3 << "****INITMINOR OTHER residual peak "<< tim.real() << LogIO::POST;
345 : //tim.mark();
346 0 : itsLoopController.setPeakResidualNoMask( peakresnomask );
347 0 : itsLoopController.setMaxPsfSidelobe( itsImages->getPSFSidelobeLevel() );
348 :
349 : //re-calculate current nsigma threhold
350 : //os<<"Calling calcRobustRMS ....syndeconv."<<LogIO::POST;
351 0 : Float nsigmathresh = 0.0;
352 0 : Bool useautomask = ( itsAutoMaskAlgorithm=="multithresh" ? true : false);
353 0 : Int iterdone = itsLoopController.getIterDone();
354 :
355 : //cerr << "INIT automask " << useautomask << " alg " << itsAutoMaskAlgorithm << " sigma " << itsNsigma << endl;
356 0 : if ( itsNsigma >0.0) {
357 0 : itsMaskHandler->setPBMaskLevel(itsPBMask);
358 0 : Array<Double> medians, robustrms;
359 : // 2 cases to use existing stats.
360 : // 1. automask has run and so the image statistics record has filled
361 : // or
362 : // 2. no automask but for the first cycle but already initial calcRMS has ran to avoid duplicate
363 : //
364 :
365 : //cerr << "useauto " << useautomask << " nfields " << itsRobustStats.nfields() << " iterdone " << iterdone << endl;
366 0 : if ((useautomask && itsRobustStats.nfields()) ||
367 0 : (!useautomask && iterdone==0 && itsRobustStats.nfields()) ) {
368 0 : os <<LogIO::DEBUG1<<"automask on: check the current stats"<<LogIO::POST;
369 : //os<< "itsRobustStats nfield="<< itsRobustStats.nfields() << LogIO::POST;;
370 0 : if (itsRobustStats.isDefined("medabsdevmed")) {
371 0 : Array<Double> mads;
372 0 : itsRobustStats.get(RecordFieldId("medabsdevmed"), mads);
373 0 : os<<LogIO::DEBUG1<<"Using robust rms from automask ="<< mads*1.4826 <<LogIO::POST;
374 0 : robustrms = mads*1.4826;
375 0 : }
376 0 : else if(itsRobustStats.isDefined("robustrms")) {
377 0 : itsRobustStats.get(RecordFieldId("robustrms"), robustrms);
378 : }
379 0 : itsRobustStats.get(RecordFieldId("median"), medians);
380 :
381 : }
382 : else { // do own stats calculation
383 0 : timer.mark();
384 :
385 0 : os<<LogIO::DEBUG1<<"Calling calcRobustRMS .. "<<LogIO::POST;
386 0 : robustrms = itsImages->calcRobustRMS(medians, itsPBMask, itsFastNoise);
387 : os<< LogIO::NORMAL << "time for calcRobustRMS: real "<< timer.real() << "s ( user " << timer.user()
388 0 : <<"s, system "<< timer.system() << "s)" << LogIO::POST;
389 : //reset itsRobustStats
390 : //cerr << "medians " << medians << " pbmask " << itsPBMask << endl;
391 : try {
392 : //os<<"current content of itsRobustStats nfields=="<<itsRobustStats.nfields()<<LogIO::POST;
393 0 : itsRobustStats.define(RecordFieldId("robustrms"), robustrms);
394 0 : itsRobustStats.define(RecordFieldId("median"), medians);
395 : }
396 0 : catch(AipsError &x) {
397 0 : throw( AipsError("Error in storing the robust image statistics") );
398 0 : }
399 :
400 : //cerr << this << " DOING robust " << itsRobustStats << endl;
401 : }
402 :
403 : /***
404 : Array<Double> robustrms =kitsImages->calcRobustRMS(medians, itsPBMask, itsFastNoise);
405 : // Before the first iteration the iteration parameters have not been
406 : // set in SIMinorCycleController
407 : // Use nsigma pass to SynthesisDeconvolver directly for now...
408 : //Float nsigma = itsLoopController.getNsigma();
409 : ***/
410 : Double minval, maxval;
411 0 : IPosition minpos, maxpos;
412 : //Double maxrobustrms = max(robustrms);
413 0 : if(robustrms.empty())
414 0 : throw(AipsError("No valid values to deconvolve"));
415 :
416 0 : minMax(minval, maxval, minpos, maxpos, robustrms);
417 :
418 : //Float nsigmathresh = nsigma * (Float)robustrms(IPosition(1,0));
419 : //nsigmathresh = itsNsigma * (Float)maxrobustrms;
420 0 : String msg="";
421 0 : if (useautomask) {
422 0 : nsigmathresh = (Float)(medians(maxpos)) + itsNsigma * (Float)maxval;
423 0 : msg+=" (nsigma*rms + median)";
424 : }
425 : else {
426 0 : nsigmathresh = itsNsigma * (Float)maxval;
427 : }
428 0 : os << "Current nsigma threshold (maximum along spectral channels ) ="<<nsigmathresh<< msg <<LogIO::POST;
429 0 : }
430 :
431 0 : itsLoopController.setNsigmaThreshold(nsigmathresh);
432 0 : itsLoopController.setPBMask(itsPBMask);
433 0 : itsLoopController.setFullSummary(itsFullSummary);
434 :
435 0 : if ( itsAutoMaskAlgorithm=="multithresh" && !initializeChanMaskFlag ) {
436 0 : IPosition maskshp = itsImages->mask()->shape();
437 0 : Int nchan = maskshp(3);
438 0 : itsChanFlag=Vector<Bool>(nchan,False);
439 0 : initializeChanMaskFlag=True;
440 : // also initialize posmask, which tracks only positive (emission)
441 :
442 0 : if(!itsPosMask){
443 : //itsPosMask = TempImage<Float> (maskshp, itsImages->mask()->coordinates(),SDMaskHandler::memoryToUse());
444 0 : itsPosMask=itsImages->tempworkimage();
445 : //you don't want to modify this here...
446 : //It is set to 0.0 in SIImageStore first time it is created.
447 : //itsPosMask->set(0);
448 0 : itsPosMask->unlock();
449 : }
450 0 : }
451 0 : os<<LogIO::DEBUG1<<"itsChanFlag.shape="<<itsChanFlag.shape()<<LogIO::POST;
452 :
453 : /*
454 : Array<Double> rmss = itsImages->calcRobustRMS();
455 : AlwaysAssert( rmss.shape()[0]>0 , AipsError);
456 : cout << "madRMS : " << rmss << " with shape : " << rmss.shape() << endl;
457 : //itsLoopController.setMadRMS( rmss[0] );
458 : */
459 :
460 0 : if( itsMaskSum != masksum || masksum == 0.0 ) // i.e. mask has changed.
461 : {
462 0 : itsMaskSum = masksum;
463 0 : itsLoopController.setMaskSum( masksum );
464 : }
465 : else // mask has not changed...
466 : {
467 0 : itsLoopController.setMaskSum( -1.0 );
468 : }
469 0 : float psfsidelobelevel = itsImages->getPSFSidelobeLevel();
470 :
471 0 : returnRecord = itsLoopController.getCycleInitializationRecord();
472 :
473 : // cerr << "INIT record " << returnRecord << endl;
474 :
475 : // itsImages->printImageStats();
476 0 : os << " Absolute Peak residual within mask : " << peakresinmask << ", over full image : " << peakresnomask << LogIO::POST;
477 0 : itsImages->releaseLocks();
478 :
479 0 : os << LogIO::DEBUG2 << "Initialized minor cycle. Returning returnRec" << LogIO::POST;
480 :
481 0 : } catch(AipsError &x) {
482 0 : throw( AipsError("Error initializing the Minor Cycle for " + itsImageName + " : "+x.getMesg()) );
483 0 : }
484 :
485 0 : return returnRecord;
486 0 : }
487 0 : void SynthesisDeconvolver::setChanFlag(const Vector<Bool>& chanflag){
488 : //ignore if it has not been given a size yet in initminorcycle
489 0 : if(itsChanFlag.nelements()==0)
490 0 : return;
491 0 : if(itsChanFlag.nelements() != chanflag.nelements())
492 0 : throw(AipsError("cannot set chan flags for different number of channels"));
493 0 : itsChanFlag =chanflag;
494 :
495 : }
496 0 : Vector<Bool> SynthesisDeconvolver::getChanFlag(){
497 0 : return itsChanFlag;
498 : }
499 0 : void SynthesisDeconvolver::setRobustStats(const Record& rec){
500 0 : itsRobustStats=Record();
501 0 : itsRobustStats=rec;
502 :
503 0 : }
504 0 : Record SynthesisDeconvolver::getRobustStats(){
505 0 : return itsRobustStats;
506 : }
507 :
508 0 : Record SynthesisDeconvolver::interactiveGUI(Record& iterRec)
509 : {
510 0 : LogIO os( LogOrigin("SynthesisDeconvolver","interactiveGUI",WHERE) );
511 0 : Record returnRecord;
512 :
513 : try {
514 :
515 : // Check that required parameters are present in the iterRec.
516 0 : Int niter=0,cycleniter=0,iterdone;
517 0 : Float threshold=0.0, cyclethreshold=0.0;
518 0 : if( iterRec.isDefined("niter") &&
519 0 : iterRec.isDefined("cycleniter") &&
520 0 : iterRec.isDefined("threshold") &&
521 0 : iterRec.isDefined("cyclethreshold") &&
522 0 : iterRec.isDefined("iterdone") )
523 : {
524 0 : iterRec.get("niter", niter);
525 0 : iterRec.get("cycleniter", cycleniter);
526 0 : iterRec.get("threshold", threshold);
527 0 : iterRec.get("cyclethreshold", cyclethreshold);
528 0 : iterRec.get("iterdone",iterdone);
529 : }
530 0 : else throw(AipsError("SD::interactiveGui() needs valid niter, cycleniter, threshold to start up."));
531 :
532 0 : if( ! itsImages ) itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
533 :
534 : // SDMaskHandler masker;
535 0 : String strthresh = String::toString(threshold)+"Jy";
536 0 : String strcycthresh = String::toString(cyclethreshold)+"Jy";
537 :
538 0 : if( itsMaskString.length()>0 ) {
539 0 : itsMaskHandler->fillMask( itsImages, itsMaskString );
540 : }
541 :
542 0 : Int iterleft = niter - iterdone;
543 0 : if( iterleft<0 ) iterleft=0;
544 :
545 0 : Int stopcode = itsMaskHandler->makeInteractiveMask( itsImages, iterleft, cycleniter, strthresh, strcycthresh );
546 :
547 0 : Quantity qa;
548 0 : casacore::Quantity::read(qa,strthresh);
549 0 : threshold = qa.getValue(Unit("Jy"));
550 0 : casacore::Quantity::read(qa,strcycthresh);
551 0 : cyclethreshold = qa.getValue(Unit("Jy"));
552 :
553 0 : itsIsMaskLoaded=true;
554 :
555 0 : returnRecord.define( RecordFieldId("actioncode"), stopcode );
556 0 : returnRecord.define( RecordFieldId("niter"), iterdone + iterleft );
557 0 : returnRecord.define( RecordFieldId("cycleniter"), cycleniter );
558 0 : returnRecord.define( RecordFieldId("threshold"), threshold );
559 0 : returnRecord.define( RecordFieldId("cyclethreshold"), cyclethreshold );
560 :
561 0 : } catch(AipsError &x) {
562 0 : throw( AipsError("Error in Interactive GUI : "+x.getMesg()) );
563 0 : }
564 0 : return returnRecord;
565 0 : }
566 :
567 :
568 0 : void SynthesisDeconvolver::setMinorCycleControl(const Record& minorCycleControlRec){
569 : //Don't know what itsloopcontroller does not need a const record;
570 0 : Record lala=minorCycleControlRec;
571 0 : itsLoopController.setCycleControls(lala);
572 :
573 0 : }
574 :
575 0 : Record SynthesisDeconvolver::executeMinorCycle(Record& minorCycleControlRec)
576 : {
577 : // LogIO os( LogOrigin("SynthesisDeconvolver","executeMinorCycle",WHERE) );
578 :
579 :
580 : // itsImages->printImageStats();
581 0 : SynthesisUtilMethods::getResource("Start Deconvolver");
582 : ///if cube execute cube deconvolution...check on residual shape as itsimagestore return 0 shape sometimes
583 0 : if(!itsImages)
584 0 : throw(AipsError("Initminor Cycle has not been called yet"));
585 0 : if(itsImages->residual()->shape()[3]> 1){
586 0 : return executeCubeMinorCycle(minorCycleControlRec);
587 : }
588 :
589 : // os << "---------------------------------------------------- Run Minor Cycle Iterations ---------------------------------------------" << LogIO::POST;
590 0 : return executeCoreMinorCycle(minorCycleControlRec);
591 : SynthesisUtilMethods::getResource("End Deconvolver");
592 : }
593 0 : Record SynthesisDeconvolver::executeCoreMinorCycle(Record& minorCycleControlRec)
594 : {
595 :
596 0 : Record returnRecord;
597 : try {
598 : // if ( !itsIsInteractive ) setAutoMask();
599 : //cerr << "MINORCYCLE control Rec " << minorCycleControlRec << endl;
600 0 : itsLoopController.setCycleControls(minorCycleControlRec);
601 0 : bool automaskon (false);
602 0 : if (itsAutoMaskAlgorithm=="multithresh") {
603 0 : automaskon=true;
604 : }
605 : //itsDeconvolver->deconvolve( itsLoopController, itsImages, itsDeconvolverId, automaskon, itsFastNoise );
606 : // include robust stats rec
607 0 : itsDeconvolver->deconvolve( itsLoopController, itsImages, itsDeconvolverId, automaskon, itsFastNoise, itsRobustStats, itsFullSummary );
608 :
609 0 : returnRecord = itsLoopController.getCycleExecutionRecord();
610 :
611 : //scatterModel(); // This is a no-op for the single-node case.
612 :
613 0 : itsImages->releaseLocks();
614 :
615 0 : } catch(AipsError &x) {
616 0 : throw( AipsError("Error in running Minor Cycle : "+x.getMesg()) );
617 0 : }
618 :
619 :
620 :
621 0 : return returnRecord;
622 0 : }
623 0 : Record SynthesisDeconvolver::executeCubeMinorCycle(Record& minorCycleControlRec, const Int AutoMaskFlag)
624 : {
625 0 : LogIO os( LogOrigin("SynthesisDeconvolver","executeCubeMinorCycle",WHERE) );
626 0 : Record returnRecord;
627 0 : Int doAutoMask=AutoMaskFlag;
628 :
629 0 : SynthesisUtilMethods::getResource("Start Deconvolver");
630 :
631 : try {
632 : // if ( !itsIsInteractive ) setAutoMask();
633 0 : if(doAutoMask < 1){
634 0 : itsLoopController.setCycleControls(minorCycleControlRec);
635 : }
636 0 : else if(doAutoMask==1){
637 0 : minorCycleControlRec=itsPreviousIterBotRec_p;
638 : }
639 0 : CubeMinorCycleAlgorithm cmc;
640 : //casa::applicator.defineAlgorithm(cmc);
641 : ///argv and argc are needed just to callthe right overloaded init
642 0 : Int argc=1;
643 0 : char **argv=nullptr;
644 0 : applicator.init(argc, argv);
645 0 : if(applicator.isController()){
646 0 : os << ((AutoMaskFlag != 1) ? "---------------------------------------------------- Run Minor Cycle Iterations ---------------------------------------------" : "---------------------------------------------------- Run Automask ---------------------------------------------" )<< LogIO::POST;
647 : /*{///TO BE REMOVED
648 : LatticeExprNode le( sum( *(itsImages->mask()) ) );
649 : os << LogIO::WARN << "#####Sum of mask BEFORE minor cycle " << le.getFloat() << endl;
650 : }
651 : */
652 0 : Timer tim;
653 0 : tim.mark();
654 : //itsImages->printImageStats();
655 : // Add itsIterdone to be sent to child processes ...needed for automask
656 : //cerr << "before record " << itsIterDone << " loopcontroller " << itsLoopController.getIterDone() << endl;
657 0 : minorCycleControlRec.define("iterdone", itsIterDone);
658 0 : if(doAutoMask < 0) // && itsPreviousIterBotRec_p.nfields() >0)
659 0 : doAutoMask=0;
660 0 : minorCycleControlRec.define("onlyautomask",doAutoMask);
661 0 : if(itsPosMask){
662 0 : minorCycleControlRec.define("posmaskname", itsPosMask->name());
663 : }
664 : //Int numprocs = applicator.numProcs();
665 : //cerr << "Number of procs: " << numprocs << endl;
666 :
667 0 : Int numchan=itsImages->residual()->shape()[3];
668 0 : Vector<Int> startchans;
669 0 : Vector<Int> endchans;
670 0 : Int numblocks=numblockchans(startchans, endchans);
671 0 : String psfname=itsImages->psf()->name();
672 :
673 0 : Float psfsidelobelevel=itsImages->getPSFSidelobeLevel();
674 0 : String residualname=itsImages->residual()->name();
675 0 : String maskname=itsImages->mask()->name();
676 0 : String modelname=itsImages->model()->name();
677 : ////Need the pb too as calcrobustrms in SynthesisDeconvolver uses it
678 0 : String pbname="";
679 0 : if(itsImages->hasPB())
680 0 : pbname=itsImages->pb()->name();
681 0 : itsImages->releaseLocks();
682 : ///in lieu of = operator go via record
683 : // need to create a proper = operator for SynthesisParamsDeconv
684 0 : SynthesisParamsDeconv decpars;
685 : ///Will have to create a = operator...right now kludging
686 : ///from record has a check that has to be bypassed for just the
687 : /// usage as a = operator
688 : {
689 0 : String tempMaskString= itsDecPars.maskString;
690 0 : itsDecPars.maskString="";
691 0 : decpars.fromRecord(itsDecPars.toRecord());
692 : //return itsDecPars back to its original state
693 0 : itsDecPars.maskString=tempMaskString;
694 0 : }
695 : ///remove starting model as already dealt with in this deconvolver
696 0 : decpars.startModel="";
697 : ///masking is dealt already by this deconvolver so mask image
698 : //is all that is needed which is sent as maskname to subdeconvolver
699 0 : decpars.maskString="";
700 0 : (decpars.maskList).resize();
701 0 : Record decParsRec = decpars.toRecord();
702 :
703 : /////Now we loop over channels and deconvolve each
704 : ///If we do want to do block of channels rather than 1 channel
705 : ///at a time chanRange can take that and the trigger to call this
706 : ///function in executeMinorCycle has to change.
707 0 : Int rank(0);
708 : Bool assigned;
709 0 : Bool allDone(false);
710 0 : Vector<Int> chanRange(2);
711 : //Record beamsetRec;
712 0 : Vector<Bool> retvals(numblocks, False);
713 0 : Vector<Bool> chanFlag(0);
714 0 : if(itsChanFlag.nelements()==0){
715 0 : itsChanFlag.resize(numchan);
716 0 : itsChanFlag.set(False);
717 : }
718 0 : Record chanflagRec;
719 0 : Int indexofretval=0;
720 0 : for (Int k=0; k < numblocks; ++k) {
721 : //os << LogIO::DEBUG1 << "deconvolving channel "<< k << LogIO::POST;
722 0 : assigned=casa::applicator.nextAvailProcess(cmc, rank);
723 : //cerr << "assigned "<< assigned << endl;
724 0 : while(!assigned) {
725 : //cerr << "SErial ? " << casa::applicator.isSerial() << endl;
726 0 : rank = casa::applicator.nextProcessDone(cmc, allDone);
727 : //cerr << "while rank " << rank << endl;
728 : //receiving output of CubeMinorCycleAlgorithm::put
729 : //#1
730 0 : Vector<Int> chanRangeProcessed;
731 0 : casa::applicator.get(chanRangeProcessed);
732 : //#2
733 :
734 0 : Record chanflagRec;
735 0 : casa::applicator.get(chanflagRec);
736 :
737 : //#3
738 0 : Record retval;
739 0 : casa::applicator.get(retval);
740 :
741 0 : Vector<Bool> retchanflag;
742 0 : chanflagRec.get("chanflag", retchanflag);
743 0 : if(retchanflag.nelements() >0)
744 0 : itsChanFlag(Slice(chanRangeProcessed[0], chanRangeProcessed[1]-chanRangeProcessed[0]+1))=retchanflag;
745 0 : Record substats=chanflagRec.asRecord("statsrec");
746 0 : setSubsetRobustStats(substats, chanRangeProcessed[0], chanRangeProcessed[1], numchan);
747 :
748 0 : retvals(indexofretval)=(retval.nfields() > 0);
749 0 : ++indexofretval;
750 : ///might need to merge these retval
751 0 : if(doAutoMask <1)
752 0 : mergeReturnRecord(retval, returnRecord, chanRangeProcessed[0]);
753 : /*if(retval.nfields())
754 : cerr << k << "deconv rank " << rank << " successful " << endl;
755 : else
756 : cerr << k << "deconv rank " << rank << " failed " << endl;
757 : */
758 : //cerr <<"rank " << rank << " return rec "<< retval << endl;
759 0 : assigned = casa::applicator.nextAvailProcess(cmc, rank);
760 :
761 0 : }
762 :
763 : ///send process info
764 : // put dec sel params #1
765 0 : applicator.put(decParsRec);
766 : // put itercontrol params #2
767 0 : applicator.put(minorCycleControlRec);
768 : // put which channel to process #3
769 0 : chanRange[0]=startchans[k]; chanRange[1]=endchans[k];
770 0 : applicator.put(chanRange);
771 : // psf #4
772 0 : applicator.put(psfname);
773 : // residual #5
774 0 : applicator.put(residualname);
775 : // model #6
776 0 : applicator.put(modelname);
777 : // mask #7
778 0 : applicator.put(maskname);
779 : //pb #8
780 0 : applicator.put(pbname);
781 : //#9 psf side lobe
782 0 : applicator.put(psfsidelobelevel);
783 : //# put chanflag
784 0 : chanFlag.resize();
785 0 : chanFlag=itsChanFlag(IPosition(1, chanRange[0]), IPosition(1, chanRange[1]));
786 :
787 0 : chanflagRec.define("chanflag", chanFlag);
788 0 : Record statrec=getSubsetRobustStats(chanRange[0], chanRange[1]);
789 0 : chanflagRec.defineRecord("statsrec", statrec);
790 0 : applicator.put(chanflagRec);
791 : /// Tell worker to process it
792 0 : applicator.apply(cmc);
793 :
794 0 : }
795 : // Wait for all outstanding processes to return
796 0 : rank = casa::applicator.nextProcessDone(cmc, allDone);
797 0 : while (!allDone) {
798 :
799 0 : Vector<Int> chanRangeProcessed;
800 0 : casa::applicator.get(chanRangeProcessed);
801 0 : Record chanflagRec;
802 0 : casa::applicator.get(chanflagRec);
803 0 : Record retval;
804 0 : casa::applicator.get(retval);
805 0 : Vector<Bool> retchanflag;
806 0 : chanflagRec.get("chanflag", retchanflag);
807 0 : if(retchanflag.nelements() >0)
808 0 : itsChanFlag(Slice(chanRangeProcessed[0], chanRangeProcessed[1]-chanRangeProcessed[0]+1))=retchanflag;
809 0 : Record substats=chanflagRec.asRecord("statsrec");
810 0 : setSubsetRobustStats(substats, chanRangeProcessed[0], chanRangeProcessed[1], numchan);
811 0 : retvals(indexofretval)=(retval.nfields() > 0);
812 0 : ++indexofretval;
813 0 : if(doAutoMask < 1)
814 0 : mergeReturnRecord(retval, returnRecord, chanRangeProcessed[0]);
815 0 : if(retval.nfields() >0)
816 : //cerr << "deconv remainder rank " << rank << " successful " << endl;
817 0 : cerr << "";
818 : else
819 0 : cerr << "deconv remainder rank " << rank << " failed " << endl;
820 :
821 0 : rank = casa::applicator.nextProcessDone(cmc, allDone);
822 0 : if(casa::applicator.isSerial())
823 0 : allDone=true;
824 :
825 0 : }
826 :
827 0 : if(anyEQ(retvals, False))
828 0 : throw(AipsError("one of more section of the cube failed in deconvolution"));
829 0 : if(doAutoMask < 1){
830 0 : itsLoopController.incrementMinorCycleCount(returnRecord.asInt("iterdone"));
831 0 : itsIterDone+=returnRecord.asInt("iterdone");
832 : }
833 0 : itsPreviousIterBotRec_p=Record();
834 0 : itsPreviousIterBotRec_p=minorCycleControlRec;
835 : /*{///TO BE REMOVED
836 : LatticeExprNode le( sum( *(itsImages->mask()) ) );
837 : os << LogIO::WARN << "#####Sum of mask AFTER minor cycle " << le.getFloat() << "loopcontroller iterdeconv " << itsLoopController.getIterDone() << endl;
838 : }*/
839 :
840 0 : }///end of if controller
841 : /////////////////////////////////////////////////
842 :
843 : //scatterModel(); // This is a no-op for the single-node case.
844 :
845 0 : itsImages->releaseLocks();
846 :
847 0 : } catch(AipsError &x) {
848 : //MPI_Abort(MPI_COMM_WORLD, 6);
849 0 : throw( AipsError("Error in running Minor Cycle : "+x.getMesg()) );
850 0 : }
851 :
852 0 : SynthesisUtilMethods::getResource("End Deconvolver");
853 :
854 0 : return returnRecord;
855 0 : }
856 :
857 : // Restore Image.
858 17 : void SynthesisDeconvolver::restore()
859 : {
860 34 : LogIO os( LogOrigin("SynthesisDeconvolver","restoreImage",WHERE) );
861 :
862 17 : if( ! itsImages )
863 : {
864 0 : itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
865 : }
866 :
867 17 : SynthesisUtilMethods::getResource("Restoration");
868 :
869 17 : itsDeconvolver->restore(itsImages);
870 17 : itsImages->releaseLocks();
871 :
872 17 : }
873 :
874 0 : void SynthesisDeconvolver::mergeReturnRecord(const Record& inRec, Record& outRec, const Int chan){
875 :
876 : ///Something has to be done about what is done in SIIterBot_state::mergeMinorCycleSummary if it is needed
877 : //int nSummaryFields = SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
878 : //int nSummaryFields = SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
879 0 : int nSummaryFields = !itsFullSummary ? 7 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
880 :
881 0 : Matrix<Double> summaryminor(nSummaryFields,0);
882 0 : if(outRec.isDefined("summaryminor"))
883 0 : summaryminor=Matrix<Double>(outRec.asArrayDouble("summaryminor"));
884 0 : Matrix<Double> subsummaryminor;
885 0 : if(inRec.isDefined("summaryminor"))
886 0 : subsummaryminor=Matrix<Double>(inRec.asArrayDouble("summaryminor"));
887 0 : if(subsummaryminor.nelements() !=0){
888 : ///The 6th element is supposed to be the subimage id
889 0 : subsummaryminor.row(5)= subsummaryminor.row(5)+(Double(chan));
890 0 : Matrix<Double> newsummary(nSummaryFields, summaryminor.shape()[1]+subsummaryminor.shape()[1]);
891 0 : Int ocol=0;
892 0 : for (Int col=0; col< summaryminor.shape()[1]; ++col, ++ocol)
893 0 : newsummary.column(ocol)=summaryminor.column(col);
894 0 : for (Int col=0; col< subsummaryminor.shape()[1]; ++col, ++ocol)
895 0 : newsummary.column(ocol)=subsummaryminor.column(col);
896 0 : summaryminor.resize(newsummary.shape());
897 0 : summaryminor=newsummary;
898 0 : }
899 0 : outRec.define("summaryminor", summaryminor);
900 : //cerr << "inRec summ minor " << inRec.asArrayDouble("summaryminor") << endl;
901 : //cerr << "outRec summ minor " << summaryminor << endl;
902 0 : outRec.define("iterdone", Int(inRec.asInt("iterdone")+ (outRec.isDefined("iterdone") ? outRec.asInt("iterdone"): Int(0))));
903 0 : outRec.define("maxcycleiterdone", outRec.isDefined("maxcycleiterdone") ? max(inRec.asInt("maxcycleiterdone"), outRec.asInt("maxcycleiterdone")) :inRec.asInt("maxcycleiterdone")) ;
904 :
905 0 : outRec.define("peakresidual", outRec.isDefined("peakresidual") ? max(inRec.asFloat("peakresidual"), outRec.asFloat("peakresidual")) :inRec.asFloat("peakresidual")) ;
906 :
907 : ///is not necessarily defined it seems
908 0 : Bool updatedmodelflag=False;
909 0 : if(inRec.isDefined("updatedmodelflag"))
910 0 : inRec.get("updatedmodelflag", updatedmodelflag);
911 0 : outRec.define("updatedmodelflag", outRec.isDefined("updatedmodelflag") ? updatedmodelflag || outRec.asBool("updatedmodelflag") : updatedmodelflag) ;
912 :
913 :
914 :
915 0 : }
916 : // get channel blocks
917 0 : Int SynthesisDeconvolver::numblockchans(Vector<Int>& startchans, Vector<Int>& endchans){
918 0 : Int nchan=itsImages->residual()->shape()[3];
919 : //roughly 8e6 pixel to deconvolve per lock/process is a minimum
920 0 : Int optchan= 8e6/(itsImages->residual()->shape()[0])/(itsImages->residual()->shape()[1]);
921 : // cerr << "OPTCHAN" << optchan << endl;
922 0 : if(optchan < 10) optchan=10;
923 0 : Int nproc= applicator.numProcs() < 2 ? 1 : applicator.numProcs()-1;
924 : /*if(nproc==1){
925 : startchans.resize(1);
926 : endchans.resize(1);
927 : startchans[0]=0;
928 : endchans[0]=nchan-1;
929 : return 1;
930 : }
931 : */
932 0 : Int blksize= nchan/nproc > optchan ? optchan : Int( std::floor(Float(nchan)/Float(nproc)));
933 0 : if(blksize< 1) blksize=1;
934 0 : Int nblk=Int(nchan/blksize);
935 0 : startchans.resize(nblk);
936 0 : endchans.resize(nblk);
937 0 : for (Int k=0; k < nblk; ++k){
938 0 : startchans[k]= k*blksize;
939 0 : endchans[k]=(k+1)*blksize-1;
940 : }
941 0 : if(endchans[nblk-1] < (nchan-1)){
942 0 : startchans.resize(nblk+1,True);
943 0 : startchans[nblk]=endchans[nblk-1]+1;
944 0 : endchans.resize(nblk+1,True);
945 0 : endchans[nblk]=nchan-1;
946 0 : ++nblk;
947 : }
948 : //cerr << "nblk " << nblk << " beg " << startchans << " end " << endchans << endl;
949 0 : return nblk;
950 : }
951 :
952 : // pbcor Image.
953 0 : void SynthesisDeconvolver::pbcor()
954 : {
955 0 : LogIO os( LogOrigin("SynthesisDeconvolver","pbcor",WHERE) );
956 :
957 0 : if( ! itsImages )
958 : {
959 0 : itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
960 : }
961 :
962 0 : itsDeconvolver->pbcor(itsImages);
963 0 : itsImages->releaseLocks();
964 :
965 0 : }
966 :
967 :
968 :
969 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
970 : //// Internal Functions start here. These are not visible to the tool layer.
971 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
972 :
973 17 : std::shared_ptr<SIImageStore> SynthesisDeconvolver::makeImageStore( String imagename, Bool noRequireSumwt )
974 : {
975 17 : std::shared_ptr<SIImageStore> imstore;
976 17 : if( itsDeconvolver->getAlgorithmName() == "mtmfs" )
977 0 : { imstore.reset( new SIImageStoreMultiTerm( imagename, itsDeconvolver->getNTaylorTerms(), true, noRequireSumwt ) ); }
978 : else
979 17 : { imstore.reset( new SIImageStore( imagename, true, noRequireSumwt ) ); }
980 :
981 17 : return imstore;
982 :
983 0 : }
984 :
985 :
986 : // #############################################
987 : // #############################################
988 : // #############################################
989 : // #############################################
990 :
991 : // Set a starting model.
992 0 : void SynthesisDeconvolver::setStartingModel()
993 : {
994 0 : LogIO os( LogOrigin("SynthesisDeconvolver","setStartingModel",WHERE) );
995 :
996 0 : if(itsAddedModel==true) {return;}
997 :
998 : try
999 : {
1000 :
1001 0 : if( itsStartingModelNames.nelements()>0 && itsImages )
1002 : {
1003 : // os << "Setting " << itsStartingModelNames << " as starting model for deconvolution " << LogIO::POST;
1004 0 : itsImages->setModelImage( itsStartingModelNames );
1005 : }
1006 :
1007 0 : itsAddedModel=true;
1008 :
1009 : }
1010 0 : catch(AipsError &x)
1011 : {
1012 0 : throw( AipsError("Error in setting starting model for deconvolution: "+x.getMesg()) );
1013 0 : }
1014 :
1015 0 : }
1016 :
1017 : // Set mask
1018 0 : Bool SynthesisDeconvolver::setupMask()
1019 : {
1020 :
1021 : ////Remembet this has to be called only after initMinorCycle
1022 0 : LogIO os( LogOrigin("SynthesisDeconvolver","setupMask",WHERE) );
1023 0 : if(!itsImages)
1024 0 : throw(AipsError("Initminor Cycle has not been called yet"));
1025 :
1026 0 : Bool maskchanged=False;
1027 : //debug
1028 0 : if( itsIsMaskLoaded==false ) {
1029 : // use mask(s)
1030 0 : if( itsMaskList[0] != "" || itsMaskType == "pb" || itsAutoMaskAlgorithm != "" ) {
1031 : // Skip automask for non-interactive mode.
1032 0 : if ( itsAutoMaskAlgorithm != "") { // && itsIsInteractive) {
1033 0 : os << "[" << itsImages->getName() << "] Setting up an auto-mask"<< ((itsPBMask>0.0)?" within PB mask limit ":"") << LogIO::POST;
1034 : ////For Cubes this is done in CubeMinorCycle
1035 : //cerr << this << "SETUP mask " << itsRobustStats << endl;
1036 0 : if(itsImages->residual()->shape()[3] ==1)
1037 0 : setAutoMask();
1038 0 : else if((itsImages->residual()->shape()[3] >1)){
1039 0 : Record dummy;
1040 0 : executeCubeMinorCycle(dummy, 1);
1041 0 : }
1042 : /***
1043 : if ( itsPBMask > 0.0 ) {
1044 : itsMaskHandler->autoMaskWithinPB( itsImages, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsPBMask);
1045 : }
1046 : else {
1047 : cerr<<"automask by automask.."<<endl;
1048 : itsMaskHandler->autoMask( itsImages, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam);
1049 : }
1050 : ***/
1051 : }
1052 :
1053 : //else if( itsMaskType=="user" && itsMaskList[0] != "" ) {
1054 0 : if( itsMaskType=="user" && itsMaskList[0] != "" ) {
1055 0 : os << "[" << itsImages->getName() << "] Setting up a mask from " << itsMaskList << ((itsPBMask>0.0)?" within PB mask limit ":"") << LogIO::POST;
1056 :
1057 0 : itsMaskHandler->fillMask( itsImages, itsMaskList);
1058 0 : if( itsPBMask > 0.0 ) {
1059 0 : itsMaskHandler->makePBMask(itsImages, itsPBMask, True);
1060 : }
1061 : }
1062 0 : else if( itsMaskType=="pb") {
1063 0 : os << "[" << itsImages->getName() << "] Setting up a PB based mask" << LogIO::POST;
1064 0 : itsMaskHandler->makePBMask(itsImages, itsPBMask);
1065 : }
1066 :
1067 0 : os << "----------------------------------------------------------------------------------------------------------------------------------------" << LogIO::POST;
1068 :
1069 : }
1070 : else {
1071 : // new im statistics creates an empty mask and need to take care of that case
1072 0 : Bool emptyMask(False);
1073 0 : if( itsImages->hasMask() )
1074 : {
1075 : // CAS-14203 - Check if mask is empty AND user didn't specify an empty mask
1076 0 : if (itsImages->getMaskSum()==0.0 && itsMaskList[0] != "") {
1077 0 : emptyMask=True;
1078 : }
1079 : }
1080 0 : if( ! itsImages->hasMask() || emptyMask ) // i.e. if there is no existing mask to re-use...
1081 : {
1082 0 : LatticeLocker lock1 (*(itsImages->mask()), FileLocker::Write);
1083 0 : if( itsIsInteractive ) itsImages->mask()->set(0.0);
1084 0 : else itsImages->mask()->set(1.0);
1085 0 : os << "[" << itsImages->getName() << "] Initializing new mask to " << (itsIsInteractive?"0.0 for interactive drawing":"1.0 for the full image") << LogIO::POST;
1086 0 : }
1087 : else {
1088 0 : os << "[" << itsImages->getName() << "] Initializing to existing mask" << LogIO::POST;
1089 : }
1090 : }
1091 :
1092 : // If anything other than automasking, don't re-make the mask here.
1093 0 : if ( itsAutoMaskAlgorithm == "" )
1094 0 : { itsIsMaskLoaded=true; }
1095 :
1096 : // Get the number of mask pixels (sum) and send to the logger.
1097 0 : Float masksum = itsImages->getMaskSum();
1098 0 : Float npix = (itsImages->getShape()).product();
1099 :
1100 : //Int npix2 = 20000*20000*16000*4;
1101 : //Float npix2f = 20000*20000*16000*4;
1102 :
1103 : //cout << " bigval : " << npix2 << " and " << npix2f << endl;
1104 :
1105 0 : os << "[" << itsImages->getName() << "] Number of pixels in the clean mask : " << masksum << " out of a total of " << npix << " pixels. [ " << 100.0 * masksum/npix << " % ]" << LogIO::POST;
1106 :
1107 0 : maskchanged=True;
1108 : }
1109 : else {
1110 : }
1111 :
1112 0 : itsImages->mask()->unlock();
1113 0 : return maskchanged;
1114 0 : }
1115 :
1116 0 : void SynthesisDeconvolver::setAutoMask()
1117 : {
1118 : //modify mask using automask otherwise no-op
1119 0 : if ( itsAutoMaskAlgorithm != "" ) {
1120 0 : itsIterDone += itsLoopController.getIterDone();
1121 :
1122 :
1123 :
1124 0 : Bool isThresholdReached = itsLoopController.isThresholdReached();
1125 : //cerr << this << " setAuto " << itsRobustStats << endl;
1126 0 : LogIO os( LogOrigin("SynthesisDeconvolver","setAutoMask",WHERE) );
1127 0 : os << "Generating AutoMask" << LogIO::POST;
1128 : //os << LogIO::WARN << "#####ItsIterDone value " << itsIterDone << endl;
1129 :
1130 : //os << "itsMinPercentChnage = " << itsMinPercentChange<< LogIO::POST;
1131 : //cerr << "SUMa of chanFlag before " << ntrue(itsChanFlag) << endl;
1132 0 : if ( itsPBMask > 0.0 ) {
1133 : //itsMaskHandler->autoMaskWithinPB( itsImages, itsPosMask, itsIterDone, itsChanFlag, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold,itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached, itsPBMask);
1134 : //pass robust stats
1135 0 : itsMaskHandler->autoMaskWithinPB( itsImages, *itsPosMask, itsIterDone, itsChanFlag, itsRobustStats, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold,itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached, itsPBMask);
1136 : }
1137 : else {
1138 : //itsMaskHandler->autoMask( itsImages, itsPosMask, itsIterDone, itsChanFlag,itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold, itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached );
1139 :
1140 : // pass robust stats
1141 0 : itsMaskHandler->autoMask( itsImages, *itsPosMask, itsIterDone, itsChanFlag, itsRobustStats, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold, itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached );
1142 : }
1143 : //cerr <<this << " SETAutoMask " << itsRobustStats << endl;
1144 : //cerr << "SUM of chanFlag AFTER " << ntrue(itsChanFlag) << endl;
1145 0 : }
1146 0 : }
1147 :
1148 : // check if restoring beam is reasonable
1149 17 : void SynthesisDeconvolver::checkRestoringBeam()
1150 : {
1151 34 : LogIO os( LogOrigin("SynthesisDeconvolver","checkRestoringBeam",WHERE) );
1152 : //check for a bad restoring beam
1153 17 : GaussianBeam beam;
1154 :
1155 17 : if( ! itsImages ) itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
1156 17 : ImageInfo psfInfo = itsImages->psf()->imageInfo();
1157 17 : if (psfInfo.hasSingleBeam()) {
1158 17 : beam = psfInfo.restoringBeam();
1159 : }
1160 0 : else if (psfInfo.hasMultipleBeams() && itsUseBeam=="common") {
1161 0 : beam = CasaImageBeamSet(psfInfo.getBeamSet()).getCommonBeam();
1162 : }
1163 17 : Double beammaj = beam.getMajor(Unit("arcsec"));
1164 17 : Double beammin = beam.getMinor(Unit("arcsec"));
1165 17 : if (std::isinf(beammaj) || std::isinf(beammin)) {
1166 0 : String msg;
1167 0 : if (itsUseBeam=="common") {
1168 0 : msg+="Bad restoring beam using the common beam (at least one of the beam axes has an infinite size) ";
1169 0 : throw(AipsError(msg));
1170 : }
1171 0 : }
1172 17 : itsImages->releaseLocks();
1173 17 : }
1174 :
1175 : // This is for interactive-clean.
1176 0 : void SynthesisDeconvolver::getCopyOfResidualAndMask( TempImage<Float> &/*residual*/,
1177 : TempImage<Float> &/*mask*/ )
1178 : {
1179 : // Actually all I think we need here are filenames JSK 12/12/12
1180 : // resize/shape and copy the residual image and mask image to these in/out variables.
1181 : // Allocate Memory here.
1182 0 : }
1183 0 : void SynthesisDeconvolver::setMask( TempImage<Float> &/*mask*/ )
1184 : {
1185 : // Here we will just pass in the new names
1186 : // Copy the input mask to the local main image mask
1187 0 : }
1188 0 : void SynthesisDeconvolver::setIterDone(const Int iterdone){
1189 : //cerr << "SETITERDONE iterdone " << iterdone << endl;
1190 : ///this get lost in initMinorCycle
1191 : //itsLoopController.incrementMinorCycleCount(iterdone);
1192 0 : itsIterDone=iterdone;
1193 :
1194 0 : }
1195 0 : void SynthesisDeconvolver::setPosMask(std::shared_ptr<ImageInterface<Float> > posmask){
1196 0 : itsPosMask=posmask;
1197 :
1198 0 : }
1199 :
1200 0 : auto key2Mat = [](const Record& rec, const String& key, const Int npol) {
1201 0 : Matrix<Double> tmp;
1202 : //cerr << "KEY2mat " << key <<" "<< rec.asArrayDouble(key).shape() << endl;
1203 0 : if(rec.asArrayDouble(key).shape().nelements()==1){
1204 0 : if(rec.asArrayDouble(key).shape()[0] != npol){
1205 0 : tmp.resize(1,rec.asArrayDouble(key).shape()[0]);
1206 0 : Vector<Double>tmpvec=rec.asArrayDouble(key);
1207 0 : tmp.row(0)=tmpvec;
1208 0 : }
1209 : else{
1210 0 : tmp.resize(rec.asArrayDouble(key).shape()[0],1);
1211 0 : Vector<Double>tmpvec=rec.asArrayDouble(key);
1212 0 : tmp.column(0)=tmpvec;
1213 0 : }
1214 :
1215 : }
1216 : else{
1217 0 : tmp=rec.asArrayDouble(key);
1218 : }
1219 0 : return tmp;
1220 0 : };
1221 :
1222 0 : Record SynthesisDeconvolver::getSubsetRobustStats(const Int chanBeg, const Int chanEnd){
1223 0 : Record outRec;
1224 : //cerr << "getSUB " << itsRobustStats << endl;
1225 0 : if(itsRobustStats.nfields()==0)
1226 0 : return outRec;
1227 0 : Matrix<Double> tmp;
1228 0 : std::vector<String> keys={"min", "max", "rms", "medabsdevmed", "med", "robustrms", "median"};
1229 0 : for (auto it = keys.begin() ; it != keys.end(); ++it){
1230 0 : if(itsRobustStats.isDefined(*it)){
1231 0 : tmp.resize();
1232 0 : tmp=key2Mat(itsRobustStats, *it, itsImages->residual()->shape()[2]);
1233 : /*
1234 : cerr << "size of " << *it << " " << itsRobustStats.asArrayDouble(*it).shape() << endl;
1235 : if(itsRobustStats.asArrayDouble(*it).shape().nelements()==1){
1236 : tmp.resize(1, itsRobustStats.asArrayDouble(*it).shape()[0]);
1237 : Vector<Double>tmpvec=itsRobustStats.asArrayDouble(*it);
1238 : tmp.row(0)=tmpvec;
1239 :
1240 : }
1241 : else
1242 : tmp=itsRobustStats.asArrayDouble(*it);
1243 : */
1244 : // cerr << std::setprecision(12) << tmp[chanBeg] << " bool " <<(tmp[chanBeg]> (DBL_MAX-(DBL_MAX*1e-15))) << endl;
1245 0 : if(tmp(0,chanBeg)> (DBL_MAX-(DBL_MAX*1e-15)))
1246 0 : return Record();
1247 : //cerr << "GETSUB blc "<< IPosition(2, 0, chanBeg)<< " trc " << IPosition(2, tmp.shape()[0]-1, chanEnd) << " shape " << tmp.shape() << endl;
1248 0 : outRec.define(*it, tmp(IPosition(2, 0, chanBeg), IPosition(2, tmp.shape()[0]-1, chanEnd)));
1249 : }
1250 : }
1251 :
1252 : // When itsImages->residual() is called, it (sometimes) creates a read-lock. Release that lock.
1253 0 : shared_ptr<ImageInterface<Float> > resimg = itsImages->residual();
1254 0 : itsImages->releaseImage( resimg );
1255 :
1256 : //cerr <<"chanbeg " << chanBeg << " chanend " << chanEnd << endl;
1257 : //cerr << "GETSUB " << outRec << endl;
1258 0 : return outRec;
1259 0 : }
1260 :
1261 0 : void SynthesisDeconvolver::setSubsetRobustStats(const Record& inrec, const Int chanBeg, const Int chanEnd, const Int numchan){
1262 0 : if(inrec.nfields()==0)
1263 0 : return ;
1264 0 : Matrix<Double> tmp;
1265 0 : std::vector<String> keys={"min", "max", "rms", "medabsdevmed", "med", "robustrms", "median"};
1266 :
1267 0 : for (auto it = keys.begin() ; it != keys.end(); ++it){
1268 0 : if(inrec.isDefined(*it)){
1269 0 : tmp.resize();
1270 0 : tmp=key2Mat(inrec, *it,itsImages->residual()->shape()[2] );
1271 0 : Matrix<Double> outvec;
1272 0 : if(itsRobustStats.isDefined(*it)){
1273 0 : outvec=key2Mat(itsRobustStats, *it, itsImages->residual()->shape()[2]);
1274 : }
1275 : else{
1276 0 : outvec.resize(itsImages->residual()->shape()[2], numchan);
1277 0 : outvec.set(DBL_MAX);
1278 : }
1279 :
1280 :
1281 0 : outvec(IPosition(2, 0, chanBeg), IPosition(2,outvec.shape()[0]-1, chanEnd))=tmp;
1282 0 : itsRobustStats.define(*it, outvec);
1283 0 : };
1284 : }
1285 :
1286 : // When itsImages->residual() is called, it (sometimes) creates a read-lock. Release that lock.
1287 0 : shared_ptr<ImageInterface<Float> > resimg = itsImages->residual();
1288 0 : itsImages->releaseImage( resimg );
1289 :
1290 : //cerr << "SETT " << itsRobustStats << endl;
1291 : //cerr << "SETT::ItsRobustStats " << Vector<Double>(itsRobustStats.asArrayDouble("min")) << endl;
1292 0 : }
1293 : } //# NAMESPACE CASA - END
1294 :
|
