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