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