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