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