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