Line data Source code
1 : //# SDMaskHandler.cc: Implementation of SDMaskHandler classes
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <stack>
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <components/ComponentModels/SkyComponent.h>
31 : #include <components/ComponentModels/ComponentList.h>
32 : #include <casacore/images/Images/ImageRegrid.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/images/Regions/RegionManager.h>
37 : #include <casacore/images/Regions/RegionHandler.h>
38 : #include <casacore/images/Regions/WCBox.h>
39 : #include <casacore/images/Regions/WCUnion.h>
40 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
41 : #include <imageanalysis/ImageAnalysis/ImageDecomposer.h>
42 : #include <imageanalysis/ImageAnalysis/ImageStatsCalculator.h>
43 : #include <imageanalysis/ImageAnalysis/Image2DConvolver.h>
44 : #include <imageanalysis/ImageAnalysis/ImageRegridder.h>
45 : #include <casacore/casa/OS/File.h>
46 : #include <casacore/lattices/LEL/LatticeExpr.h>
47 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
48 : #include <casacore/lattices/Lattices/LatticeStepper.h>
49 : #include <casacore/lattices/Lattices/LatticeIterator.h>
50 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
51 : #include <casacore/lattices/LRegions/LCEllipsoid.h>
52 : #include <casacore/lattices/LRegions/LCUnion.h>
53 : #include <casacore/lattices/LRegions/LCExtension.h>
54 : #include <casacore/lattices/LRegions/LCPagedMask.h>
55 : #include <synthesis/TransformMachines/StokesImageUtil.h>
56 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
57 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
58 : #include <casacore/casa/Exceptions/Error.h>
59 : #include <casacore/casa/BasicSL/String.h>
60 : #include <casacore/casa/Utilities/Assert.h>
61 : #include <casacore/casa/OS/Directory.h>
62 : #include <casacore/tables/Tables/TableLock.h>
63 :
64 : #include <sstream>
65 :
66 : #include <casacore/casa/Logging/LogMessage.h>
67 : #include <casacore/casa/Logging/LogIO.h>
68 : #include <casacore/casa/Logging/LogSink.h>
69 :
70 : #include <imageanalysis/Annotations/RegionTextList.h>
71 : #include <synthesis/ImagerObjects/SDMaskHandler.h>
72 :
73 :
74 : using namespace casacore;
75 : namespace casa { //# NAMESPACE CASA - BEGIN
76 :
77 2094 : SDMaskHandler::SDMaskHandler()
78 : {
79 2094 : itsMax = DBL_MAX;
80 2094 : itsRms = DBL_MAX;
81 2094 : itsSidelobeLevel = 0.0;
82 : //itsPBMaskLevel = 0.0;
83 : //itsTooLongForFname=false;
84 : #ifdef NAME_MAX
85 2094 : itsNAME_MAX = NAME_MAX;
86 : #else
87 : itsNAME_MAX = 255;
88 : #endif
89 :
90 2094 : itsTooLongForFname = false;
91 :
92 2094 : }
93 :
94 :
95 2094 : SDMaskHandler::~SDMaskHandler()
96 : {
97 2094 : }
98 :
99 0 : void SDMaskHandler::resetMask(std::shared_ptr<SIImageStore> imstore)
100 : {
101 0 : imstore->mask()->set(1.0);
102 0 : imstore->mask()->unlock();
103 0 : }
104 :
105 : /***
106 : void SDMaskHandler::fillMask(std::shared_ptr<SIImageStore> imstore, Vector<String> maskStrings)
107 : {
108 : String maskString;
109 : if (maskStrings.nelements()) {
110 : for (uInt imsk = 0; imsk < maskStrings.nelements(); imsk++) {
111 : maskString = maskStrings[imsk];
112 : fillMask(imstore, maskString);
113 : }
114 : }
115 : else {
116 : fillMask(imstore, String(""));
117 : }
118 : }
119 : ***/
120 :
121 62 : void SDMaskHandler::fillMask(std::shared_ptr<SIImageStore> imstore, Vector<String> maskStrings)
122 : {
123 186 : LogIO os( LogOrigin("SDMaskHandler","fillMask",WHERE) );
124 124 : String maskString;
125 : try {
126 186 : TempImage<Float> tempAllMaskImage(imstore->mask()->shape(), imstore->mask()->coordinates(), memoryToUse());
127 62 : tempAllMaskImage.set(0.0);
128 62 : if (maskStrings.nelements()) {
129 : //working temp mask image
130 124 : TempImage<Float> tempMaskImage(imstore->mask()->shape(), imstore->mask()->coordinates(), memoryToUse());
131 62 : copyMask(*(imstore->mask()), tempMaskImage);
132 124 : for (uInt imsk = 0; imsk < maskStrings.nelements(); imsk++) {
133 62 : maskString = maskStrings[imsk];
134 62 : bool isTable(false);
135 62 : bool isCasaImage(false);
136 : // test
137 124 : Path testPath(maskString);
138 62 : if (testPath.baseName()==maskString) {
139 62 : if (int(maskString.length()) > itsNAME_MAX)
140 : // the string is too long if it is file name or could be region text
141 1 : itsTooLongForFname=true;
142 : }
143 62 : if (maskString!="") {
144 62 : if(!itsTooLongForFname)
145 61 : isTable = Table::isReadable(maskString);
146 62 : if (!isTable) {
147 : try {
148 46 : if (ImageOpener::imageType(maskString)==ImageOpener::IMAGECONCAT) isCasaImage=true;
149 : //os <<"IT'S A CONCATImage"<<LogIO::POST;
150 : }
151 0 : catch (AipsError &x) {
152 0 : os <<"cannot find imageType... "<<LogIO::POST;
153 : }
154 : }
155 : else {
156 32 : Table imtab = Table(maskString, Table::Old);
157 32 : Vector<String> colnames = imtab.tableDesc().columnNames();
158 16 : if ( colnames[0]=="map" ) isCasaImage=true;
159 : }
160 62 : if (isCasaImage) {
161 16 : std::shared_ptr<ImageInterface<Float> > inmaskptr;
162 16 : LatticeBase* latt =ImageOpener::openImage(maskString);
163 16 : inmaskptr.reset(dynamic_cast<ImageInterface<Float>* >(latt));
164 32 : IPosition inShape = inmaskptr->shape();
165 32 : IPosition outShape = imstore->mask()->shape();
166 16 : Int specAxis = CoordinateUtil::findSpectralAxis(inmaskptr->coordinates());
167 16 : Int inNchan = (specAxis==-1? 1: inShape(specAxis) );
168 16 : Int outSpecAxis = CoordinateUtil::findSpectralAxis(imstore->mask()->coordinates());
169 32 : Vector <Stokes::StokesTypes> inWhichPols, outWhichPols;
170 16 : Int stokesAxis = CoordinateUtil::findStokesAxis(inWhichPols, inmaskptr->coordinates());
171 16 : Int inNstokes = (stokesAxis==-1? 1: inShape(stokesAxis) );
172 16 : Int outStokesAxis = CoordinateUtil::findStokesAxis(outWhichPols, imstore->mask()->coordinates());
173 16 : if (inNchan == 1 && outShape(outSpecAxis)>1) {
174 6 : os << "Extending mask image: " << maskString << LogIO::POST;
175 6 : expandMask(*inmaskptr, tempMaskImage);
176 : }
177 10 : else if(inNstokes == 1 && outShape(outStokesAxis) > 1 ) {
178 4 : os << "Extending mask image along Stokes axis: " << maskString << LogIO::POST;
179 4 : expandMask(*inmaskptr, tempMaskImage);
180 : }
181 : else {
182 6 : os << "Copying mask image: " << maskString << LogIO::POST;
183 6 : copyMask(*inmaskptr, tempMaskImage);
184 : }
185 : }// end of readable table
186 : else {
187 : //
188 46 : std::unique_ptr<Record> myrec(nullptr);
189 : try {
190 46 : if (!itsTooLongForFname) {
191 135 : myrec.reset(RegionManager::readImageFile(maskString,String("temprgrec")));
192 0 : if (myrec!=0) {
193 0 : Bool ret(false);
194 0 : Matrix<Quantity> dummyqmat;
195 0 : Matrix<Float> dummyfmat;
196 0 : ret=SDMaskHandler::regionToImageMask(tempMaskImage, myrec.get(), dummyqmat, dummyfmat);
197 0 : if (!ret) cout<<"regionToImageMask failed..."<<endl;
198 0 : os << "Reading region record mask: " << maskString << LogIO::POST;
199 : }
200 : }
201 : else {
202 : // skip to check the mask as a direct region text input
203 1 : throw(AipsError("Too long for filename. Try the string as a direct region text"));
204 : }
205 : }
206 92 : catch (AipsError &x) {
207 : try {
208 46 : ImageRegion* imageRegion=0;
209 46 : os << LogIO::NORMAL3<< "Reading text mask: " << maskString << LogIO::POST;
210 46 : SDMaskHandler::regionTextToImageRegion(maskString, tempMaskImage, imageRegion);
211 46 : if (imageRegion!=0)
212 46 : SDMaskHandler::regionToMask(tempMaskImage,*imageRegion, Float(1.0));
213 : }
214 0 : catch (AipsError &x) {
215 0 : os << LogIO::WARN << maskString << "is invalid mask. Skipping this mask..." << LogIO::POST;
216 : }
217 : }
218 : }// end of region string
219 : }// end of non-emtpy maskstring
220 :
221 124 : LatticeExpr<Float> addedmask( tempMaskImage + tempAllMaskImage );
222 62 : tempAllMaskImage.copyData( LatticeExpr<Float>( iif(addedmask > 0.0, 1.0, 0.0) ) );
223 : }
224 62 : imstore->mask()->copyData(tempAllMaskImage);
225 62 : imstore->mask()->unlock();
226 : }
227 0 : } catch( AipsError &x )
228 : {
229 0 : throw(AipsError("Error in constructing "+ imstore->getName() +".mask from " + maskString + " : " + x.getMesg()));
230 : }
231 62 : }
232 :
233 :
234 :
235 0 : void SDMaskHandler::fillMask(std::shared_ptr<SIImageStore> imstore, String maskString)
236 : {
237 :
238 : try {
239 :
240 : //// imstore->mask() will return a pointer to an ImageInterface (allocation happens on first access).
241 :
242 : // cout << "Call makeMask here to fill in " << imstore->mask()->name() << " from " << maskString << ". For now, set mask to 1 inside a central box" << endl;
243 :
244 : //interpret maskString
245 0 : if (maskString !="") {
246 0 : bool isTable(false);
247 0 : Path testPath(maskString);
248 0 : if (testPath.baseName() == maskString) {
249 0 : if (int(maskString.length()) > itsNAME_MAX)
250 0 : itsTooLongForFname=true;
251 : }
252 0 : if (!itsTooLongForFname)
253 0 : isTable = Table::isReadable(maskString);
254 : //if (!itsTooLongForFname && Table::isReadable(maskString) ) {
255 0 : if ( isTable ) {
256 0 : Table imtab = Table(maskString, Table::Old);
257 0 : Vector<String> colnames = imtab.tableDesc().columnNames();
258 0 : if ( colnames[0]=="map" ) {
259 :
260 : // looks like a CASA image ... probably should check coord exists in the keyword also...
261 : // cout << "copy this input mask...."<<endl;
262 0 : PagedImage<Float> inmask(maskString);
263 0 : IPosition inShape = inmask.shape();
264 0 : IPosition outShape = imstore->mask()->shape();
265 0 : Int specAxis = CoordinateUtil::findSpectralAxis(inmask.coordinates());
266 0 : Int outSpecAxis = CoordinateUtil::findSpectralAxis(imstore->mask()->coordinates());
267 0 : if (inShape(specAxis) == 1 && outShape(outSpecAxis)>1) {
268 0 : expandMask(inmask, *(imstore->mask()));
269 : }
270 : else {
271 0 : copyMask(inmask, *(imstore->mask()));
272 : }
273 : }// end of ''map''
274 : }// end of readable table
275 : else {
276 0 : ImageRegion* imageRegion=0;
277 0 : SDMaskHandler::regionTextToImageRegion(maskString, *(imstore->mask()), imageRegion);
278 0 : if (imageRegion!=0)
279 0 : SDMaskHandler::regionToMask(*(imstore->mask()),*imageRegion, Float(1.0));
280 : }// end of region string
281 : }
282 : else {
283 : /////// Temporary code to set a mask in the inner quarter.
284 : /////// This is only for testing... should go when 'maskString' can be used to fill it in properly.
285 0 : IPosition imshp = imstore->mask()->shape();
286 0 : AlwaysAssert( imshp.nelements() >=2 , AipsError );
287 :
288 0 : Slicer themask;
289 0 : IPosition blc(imshp.nelements(), 0);
290 0 : IPosition trc = imshp-1;
291 0 : IPosition inc(imshp.nelements(), 1);
292 :
293 0 : blc(0)=int(imshp[0]*0.25);
294 0 : blc(1)=int(imshp[1]*0.25);
295 0 : trc(0)=int(imshp[0]*0.75);
296 0 : trc(1)=int(imshp[1]*0.75);
297 :
298 0 : LCBox::verify(blc, trc, inc, imshp);
299 0 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
300 :
301 0 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(*(imstore->mask()), imslice, true) );
302 0 : referenceImage->set(1.0);
303 : }
304 :
305 0 : imstore->mask()->unlock();
306 :
307 0 : } catch( AipsError &x )
308 : {
309 0 : throw(AipsError("Error in constructing "+ imstore->getName() +".mask from " + maskString + " : " + x.getMesg()));
310 : }
311 0 : }
312 :
313 : //void SDMaskHandler::makeMask()
314 0 : std::shared_ptr<ImageInterface<Float> > SDMaskHandler::makeMask(const String& maskName, const Quantity threshold,
315 : //void SDMaskHandler::makeMask(const String& maskName, const Quantity threshold,
316 : ImageInterface<Float>& tempim)
317 : // ImageInterface<Float>& tempim,
318 : // ImageInterface<Float> *newMaskImage)
319 : {
320 0 : LogIO os( LogOrigin("SDMaskHandler","makeMask",WHERE) );
321 : //
322 : // create mask from a threshold... Imager::mask()...
323 : //default handling?
324 0 : String maskFileName(maskName);
325 0 : if ( maskFileName=="" ) {
326 0 : maskFileName = tempim.name() + ".mask";
327 : }
328 0 : if (!cloneImShape(tempim, maskFileName)) {
329 0 : throw(AipsError("Cannot make a mask from "+tempim.name()));
330 : }
331 0 : PagedImage<Float> *newMaskImage = new PagedImage<Float>(maskFileName, TableLock::AutoNoReadLocking);
332 : //newMaskImage = PagedImage<Float>(maskFileName, TableLock::AutoNoReadLocking);
333 : //PagedImage<Float>(maskFileName, TableLock::AutoNoReadLocking);
334 0 : StokesImageUtil::MaskFrom(*newMaskImage, tempim, threshold);
335 0 : return std::shared_ptr<ImageInterface<Float> >(newMaskImage);
336 : }
337 :
338 : //Bool SDMaskHandler::regionToImageMask(const String& maskName, Record* regionRec, Matrix<Quantity>& blctrcs,
339 0 : Bool SDMaskHandler::regionToImageMask(ImageInterface<Float>& maskImage, const Record* regionRec, const Matrix<Quantity>& blctrcs,
340 : const Matrix<Float>& circles, const Float& value) {
341 :
342 0 : LogIO os(LogOrigin("imager", "regionToImageMask", WHERE));
343 :
344 : try {
345 : //PagedImage<Float> tempmask(TiledShape(maskImage->shape(),
346 : // maskImage->niceCursorShape()), maskImage->coordinates(), tempfname);
347 0 : std::shared_ptr<ImageInterface<Float> > tempmask;
348 0 : tempmask.reset( new TempImage<Float>(TiledShape(maskImage.shape(),maskImage.niceCursorShape()), maskImage.coordinates(), memoryToUse() ) );
349 : //tempmask = new PagedImage<Float>(maskImage.shape(), maskImage.coordinates(),"__tmp_rgmask");
350 0 : tempmask->copyData(maskImage);
351 :
352 0 : CoordinateSystem cSys=tempmask->coordinates();
353 : //maskImage.table().markForDelete();
354 0 : ImageRegion *boxregions=0;
355 0 : ImageRegion *circleregions=0;
356 0 : RegionManager regMan;
357 0 : regMan.setcoordsys(cSys);
358 0 : if (blctrcs.nelements()!=0){
359 0 : boxRegionToImageRegion(*tempmask, blctrcs, boxregions);
360 : }
361 0 : if (circles.nelements()!=0) {
362 0 : circleRegionToImageRegion(*tempmask, circles, circleregions);
363 : }
364 0 : ImageRegion* recordRegion=0;
365 0 : if(regionRec !=0){
366 : //if(regionRec.nfields() !=0){
367 0 : recordRegionToImageRegion(regionRec, recordRegion);
368 : }
369 :
370 0 : ImageRegion *unionReg=0;
371 0 : if(boxregions!=0 && recordRegion!=0){
372 0 : unionReg=regMan.doUnion(*boxregions, *recordRegion);
373 0 : delete boxregions; boxregions=0;
374 0 : delete recordRegion; recordRegion=0;
375 : }
376 0 : else if(boxregions !=0){
377 0 : unionReg=boxregions;
378 : }
379 0 : else if(recordRegion !=0){
380 0 : unionReg=recordRegion;
381 : }
382 :
383 0 : if(unionReg !=0){
384 0 : regionToMask(*tempmask, *unionReg, value);
385 0 : delete unionReg; unionReg=0;
386 : }
387 : //As i can't unionize LCRegions and WCRegions; do circles seperately
388 0 : if(circleregions !=0){
389 0 : regionToMask(*tempmask, *circleregions, value);
390 0 : delete circleregions;
391 0 : circleregions=0;
392 : }
393 : //maskImage.table().unmarkForDelete();
394 0 : maskImage.copyData(*tempmask);
395 : }
396 0 : catch (AipsError& x) {
397 0 : os << "Error in regionToMaskImage() : " << x.getMesg() << LogIO::EXCEPTION;
398 : }
399 0 : return true;
400 : }
401 :
402 46 : Bool SDMaskHandler::regionToMask(ImageInterface<Float>& maskImage, ImageRegion& imageregion, const Float& value)
403 : {
404 92 : SubImage<Float> partToMask(maskImage, imageregion, true);
405 92 : LatticeRegion latReg=imageregion.toLatticeRegion(maskImage.coordinates(), maskImage.shape());
406 92 : ArrayLattice<Bool> pixmask(latReg.get());
407 92 : LatticeExpr<Float> myexpr(iif(pixmask, value, partToMask) );
408 46 : partToMask.copyData(myexpr);
409 92 : return true;
410 : }
411 :
412 0 : void SDMaskHandler::boxRegionToImageRegion(const ImageInterface<Float>& maskImage, const Matrix<Quantity>& blctrcs, ImageRegion*& boxImageRegions)
413 : {
414 0 : if(blctrcs.shape()(1) != 4)
415 0 : throw(AipsError("Need a list of 4 elements to define a box"));
416 :
417 0 : CoordinateSystem cSys=maskImage.coordinates();
418 0 : RegionManager regMan;
419 0 : regMan.setcoordsys(cSys);
420 0 : Vector<Quantum<Double> > blc(2);
421 0 : Vector<Quantum<Double> > trc(2);
422 0 : Int nrow=blctrcs.shape()(0);
423 0 : Vector<Int> absRel(2, RegionType::Abs);
424 0 : PtrBlock<const WCRegion *> lesbox(nrow);
425 0 : for (Int k=0; k < nrow; ++k){
426 0 : blc(0) = blctrcs(k,0);
427 0 : blc(1) = blctrcs(k,1);
428 0 : trc(0) = blctrcs(k,2);
429 0 : trc(1) = blctrcs(k,3);
430 0 : lesbox[k]= new WCBox (blc, trc, cSys, absRel);
431 : }
432 0 : boxImageRegions=regMan.doUnion(lesbox);
433 0 : if (boxImageRegions!=0) {
434 : }
435 0 : for (Int k=0; k < nrow; ++k){
436 0 : delete lesbox[k];
437 : }
438 0 : }
439 :
440 0 : void SDMaskHandler::circleRegionToImageRegion(const ImageInterface<Float>& maskImage, const Matrix<Float>& circles,
441 : ImageRegion*& circleImageRegions)
442 : {
443 0 : if(circles.shape()(1) != 3)
444 0 : throw(AipsError("Need a list of 3 elements to define a circle"));
445 :
446 0 : CoordinateSystem cSys=maskImage.coordinates();
447 0 : RegionManager regMan;
448 0 : regMan.setcoordsys(cSys);
449 0 : Int nrow=circles.shape()(0);
450 0 : Vector<Float> cent(2);
451 0 : cent(0)=circles(0,1); cent(1)=circles(0,2);
452 0 : Float radius=circles(0,0);
453 0 : IPosition xyshape(2,maskImage.shape()(0),maskImage.shape()(1));
454 0 : LCEllipsoid *circ= new LCEllipsoid(cent, radius, xyshape);
455 : //Tell LCUnion to delete the pointers
456 0 : LCUnion *elunion= new LCUnion(true, circ);
457 : //now lets do the remainder
458 0 : for (Int k=1; k < nrow; ++k){
459 0 : cent(0)=circles(k,1); cent(1)=circles(k,2);
460 0 : radius=circles(k,0);
461 0 : circ= new LCEllipsoid(cent, radius, xyshape);
462 0 : elunion=new LCUnion(true, elunion, circ);
463 : }
464 : //now lets extend that to the whole image
465 0 : IPosition trc(2);
466 0 : trc(0)=maskImage.shape()(2)-1;
467 0 : trc(1)=maskImage.shape()(3)-1;
468 0 : LCBox lbox(IPosition(2,0,0), trc,
469 0 : IPosition(2,maskImage.shape()(2),maskImage.shape()(3)) );
470 0 : LCExtension linter(*elunion, IPosition(2,2,3),lbox);
471 0 : circleImageRegions=new ImageRegion(linter);
472 0 : delete elunion;
473 0 : }
474 :
475 46 : void SDMaskHandler::recordRegionToImageRegion(const Record* imageRegRec, ImageRegion*& imageRegion )
476 : //void SDMaskHandler::recordRegionToImageRegion(Record& imageRegRec, ImageRegion*& imageRegion )
477 : {
478 46 : if(imageRegRec !=0){
479 46 : TableRecord rec1;
480 46 : rec1.assign(*imageRegRec);
481 46 : imageRegion=ImageRegion::fromRecord(rec1,"");
482 : }
483 46 : }
484 :
485 :
486 46 : void SDMaskHandler::regionTextToImageRegion(const String& text, const ImageInterface<Float>& regionImage,
487 : ImageRegion*& imageRegion)
488 : {
489 138 : LogIO os( LogOrigin("SDMaskHandler", "regionTextToImageRegion",WHERE) );
490 :
491 : try {
492 92 : IPosition imshape = regionImage.shape();
493 92 : CoordinateSystem csys = regionImage.coordinates();
494 92 : File fname;
495 : int maxFileName;
496 46 : bool skipfnamecheck(false);
497 : #ifdef NAME_MAX
498 46 : maxFileName=NAME_MAX;
499 : #else
500 : maxFileName=255;
501 : #endif
502 :
503 92 : Path testPath(text);
504 46 : if (testPath.baseName() == text) {
505 46 : if (int(text.length()) > maxFileName)
506 : // could be direct region text
507 1 : skipfnamecheck=true;
508 : }
509 46 : fname=text;
510 46 : Record* imageRegRec=0;
511 92 : Record myrec;
512 :
513 : //Record imageRegRec;
514 46 : if (!skipfnamecheck && fname.exists() && fname.isRegular()) {
515 50 : RegionTextList CRTFList(text, csys, imshape);
516 25 : myrec = CRTFList.regionAsRecord();
517 : }
518 : else { // direct text input....
519 63 : Regex rx (Regex::fromPattern("*\\[*\\]*"));
520 21 : if (text.matches(rx)) {
521 42 : RegionTextList CRTFList(csys, text, imshape);
522 21 : myrec = CRTFList.regionAsRecord();
523 : //cerr<<"myrec.nfields()="<<myrec.nfields()<<endl;
524 : }
525 : else {
526 0 : throw(AipsError("Input mask, '"+text+"' does not exist if it is inteded as a mask file name."+
527 0 : "Or invalid CRTF syntax if it is intended as a direct region specification."));
528 : }
529 : }
530 46 : imageRegRec = new Record();
531 46 : imageRegRec->assign(myrec);
532 46 : recordRegionToImageRegion(imageRegRec, imageRegion);
533 46 : delete imageRegRec;
534 : }
535 0 : catch (AipsError& x) {
536 0 : os << LogIO::SEVERE << "Exception: "<< x.getMesg() << LogIO::POST;
537 : }
538 46 : }
539 :
540 0 : void SDMaskHandler::copyAllMasks(const Vector< std::shared_ptr<ImageInterface<Float> > > inImageMasks, ImageInterface<Float>& outImageMask)
541 : {
542 0 : LogIO os( LogOrigin("SDMaskHandler", "copyAllMasks", WHERE) );
543 :
544 0 : TempImage<Float> tempoutmask(outImageMask.shape(), outImageMask.coordinates(), memoryToUse());
545 :
546 0 : for (uInt i = 0; i < inImageMasks.nelements(); i++) {
547 0 : copyMask(*inImageMasks(i), tempoutmask);
548 0 : outImageMask.copyData( (LatticeExpr<Float>)(tempoutmask+outImageMask) );
549 : }
550 0 : }
551 :
552 68 : void SDMaskHandler::copyMask(const ImageInterface<Float>& inImageMask, ImageInterface<Float>& outImageMask)
553 : {
554 204 : LogIO os( LogOrigin("SDMaskHandler", "copyMask", WHERE) );
555 :
556 : //output mask coords
557 136 : IPosition outshape = outImageMask.shape();
558 136 : CoordinateSystem outcsys = outImageMask.coordinates();
559 136 : DirectionCoordinate outDirCsys = outcsys.directionCoordinate();
560 136 : SpectralCoordinate outSpecCsys = outcsys.spectralCoordinate();
561 :
562 : // do regrid
563 136 : IPosition axes(3,0,1,2);
564 136 : IPosition inshape = inImageMask.shape();
565 136 : CoordinateSystem incsys = inImageMask.coordinates();
566 136 : DirectionCoordinate inDirCsys = incsys.directionCoordinate();
567 136 : SpectralCoordinate inSpecCsys = incsys.spectralCoordinate();
568 : //Check the conversion layer consistentcy between input and output.
569 : //Change the frame of the convesion layer in incsys to that of outcsys if different.
570 68 : if (inSpecCsys.frequencySystem(True)!=outSpecCsys.frequencySystem(True)) {
571 0 : incsys.setSpectralConversion(MFrequency::showType(outSpecCsys.frequencySystem(True)));
572 : }
573 :
574 136 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(incsys);
575 68 : axes(0) = dirAxes(0);
576 68 : axes(1) = dirAxes(1);
577 68 : axes(2) = CoordinateUtil::findSpectralAxis(incsys);
578 :
579 : //const String outfilename = outImageMask.name()+"_"+String::toString(HostInfo::processID());
580 :
581 : try {
582 : // Since regrid along the spectral axis does not seem to work
583 : // properly, replacing with ImageRegridder
584 : //ImageRegrid<Float> imregrid;
585 : //imregrid.showDebugInfo(1);
586 : //imregrid.regrid(outImageMask, Interpolate2D::LINEAR, axes, inImageMask);
587 : //
588 68 : TempImage<Float>* inImageMaskptr = new TempImage<Float>(inshape,incsys,memoryToUse());
589 136 : ArrayLattice<Bool> inimmask(inImageMask.getMask());
590 68 : inImageMaskptr->copyData((LatticeExpr<Float>)(inImageMask * iif(inimmask,1.0,0.0)) );
591 : //
592 136 : SPCIIF tempim(inImageMaskptr);
593 136 : SPCIIF templateim(new TempImage<Float>(outshape,outcsys, memoryToUse()));
594 68 : Record* dummyrec = 0;
595 : //ImageRegridder regridder(tempim, outfilename, templateim, axes, dummyrec, "", true, outshape);
596 : //TTDEBUG
597 204 : ImageRegridder<Float> regridder(tempim, "", templateim, axes, dummyrec, "", true, outshape);
598 68 : regridder.setMethod(Interpolate2D::LINEAR);
599 136 : SPIIF retim = regridder.regrid();
600 : //outImageMask.copyData( (LatticeExpr<Float>) iif(*retim > 0.1, 1.0, 0.0) );
601 68 : ArrayLattice<Bool> retimmask(retim->getMask());
602 : //LatticeExpr<Float> withmask( (*retim) * iif(retimmask,1.0,0.0) );
603 : //outImageMask.copyData( withmask );
604 68 : outImageMask.copyData( (LatticeExpr<Float>)((*retim) * iif(retimmask,1.0,0.0)) );
605 : //LatticeUtilities::copyDataAndMask(os, outImageMask, *retim );
606 : //} catch (AipsError &x) {
607 : // throw(AipsError("Image regrid error : "+ x.getMesg()));
608 : // }
609 0 : } catch (AipsError &x) {
610 0 : os<<LogIO::WARN <<x.getMesg()<<LogIO::POST;
611 0 : os<<LogIO::WARN <<" Cannot regrid the input mask to output mask, it will be empty."<<LogIO::POST;
612 0 : outImageMask.set(0);
613 : }
614 :
615 : // no longer needed (output file in regrid is now set to "" so no need of this clean-up)
616 : //try
617 : // {
618 : // delete the outfilename image on disk
619 : // Directory dd(outfilename);
620 : // dd.removeRecursive();
621 : // }
622 : //catch (AipsError &x) {
623 : // throw(AipsError("Cannot delete temporary mask image : "+ x.getMesg()));
624 : // os << LogIO::WARN << "Cannot delete temporary mask image : " << x.getMesg() << LogIO::POST;
625 : //}
626 :
627 68 : }
628 :
629 10 : void SDMaskHandler::expandMask(const ImageInterface<Float>& inImageMask, ImageInterface<Float>& outImageMask)
630 : {
631 30 : LogIO os( LogOrigin("SDMaskHandler", "expandMask", WHERE) );
632 :
633 : // expand mask with input range (in spectral axis and stokes?) ... to output range on outimage
634 : // currently expand a continuum mask to a cube mask in channels (to all channels)
635 : // or continuum Stokes I mask to multi-Stokes mask
636 : // or continuum with multi-Stokes mask to cube with multi-Stokes mask
637 20 : IPosition inShape = inImageMask.shape();
638 20 : CoordinateSystem inCsys = inImageMask.coordinates();
639 20 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(inCsys);
640 10 : Int inSpecAxis = CoordinateUtil::findSpectralAxis(inCsys);
641 : Int inNchan;
642 10 : if (inSpecAxis==-1) {
643 4 : inNchan=1;
644 : }
645 : else {
646 6 : inNchan = inShape(inSpecAxis);
647 : }
648 :
649 20 : Vector<Stokes::StokesTypes> inWhichPols;
650 10 : Int inStokesAxis = CoordinateUtil::findStokesAxis(inWhichPols,inCsys);
651 : Int inNpol;
652 10 : if (inStokesAxis==-1) {
653 4 : inNpol=1;
654 : }
655 : else {
656 6 : inNpol = inShape(inStokesAxis);
657 : }
658 :
659 : //
660 20 : IPosition outShape = outImageMask.shape();
661 20 : CoordinateSystem outCsys = outImageMask.coordinates();
662 20 : Vector<Int> outDirAxes = CoordinateUtil::findDirectionAxes(outCsys);
663 10 : Int outSpecAxis = CoordinateUtil::findSpectralAxis(outCsys);
664 10 : Int outNchan = outShape(outSpecAxis);
665 20 : Vector<Stokes::StokesTypes> outWhichPols;
666 10 : Int outStokesAxis = CoordinateUtil::findStokesAxis(outWhichPols,outCsys);
667 10 : Int outNpol = outShape(outStokesAxis);
668 :
669 20 : IPosition start(4,0,0,0,0);
670 20 : IPosition length(4,outShape(outDirAxes(0)), outShape(outDirAxes(1)),1,1);
671 : //Slicer sl(start, length);
672 :
673 10 : Int stokesInc = 1;
674 : // for fixing removed degenerate axis
675 10 : Bool addSpecAxis = (inSpecAxis == -1? True: False);
676 : // Do expansion for input mask with single channel (continuum)
677 10 : if (inNchan==1 ) {
678 8 : if (inNpol == 1 ) {
679 6 : stokesInc = 1;
680 : }
681 2 : else if (inShape(inStokesAxis)==outShape(outStokesAxis)) {
682 2 : stokesInc = inShape(inStokesAxis);
683 : }
684 : else {
685 0 : throw(AipsError("Cannot extend the input mask of "+String::toString(inNpol)+
686 0 : " Stokes planes to the output mask of "+String::toString(outNpol)+
687 0 : " Stokes planes. Please modify the input mask to make it a Stokes I mask or a mask with the same number of Stokes planes as the output.") );
688 : }
689 :
690 8 : length(outStokesAxis) = stokesInc;
691 :
692 : // I stokes cont -> cube: regrid ra.dec on the input single plane
693 : // I stokes cont -> cont multi-stokes: regrid ra.dec on the input
694 : // I stokes cont -> cube multi-stokes: regid ra.dec on input
695 :
696 16 : Slicer sl(start, length);
697 : // make a subImage for regridding output
698 16 : SubImage<Float> chanMask(outImageMask, sl, true);
699 16 : ImageRegrid<Float> imregrid;
700 24 : TempImage<Float> tempInImageMask(chanMask.shape(), chanMask.coordinates(),memoryToUse());
701 8 : std::unique_ptr<ImageInterface<Float> > outmaskim_ptr;
702 8 : if ( chanMask.shape().nelements() > inImageMask.shape().nelements() ) {
703 8 : String stokesStr;
704 4 : if (inNpol==1) {
705 3 : stokesStr="I";
706 : }
707 : else {
708 1 : stokesStr="";
709 : //for (uInt ipol=0; ipol < inWhichPols.nelements(); ipol++) {
710 : // stokesStr+=Stokes::name(inWhichPols(ipol));
711 : //}
712 : }
713 : //os<<"Adding degenerate axes: addSpecAxis="<<addSpecAxis<<" stokes="<<stokesStr<<LogIO::POST;
714 4 : casacore::ImageUtilities::addDegenerateAxes(os, outmaskim_ptr, inImageMask,"",False, addSpecAxis, stokesStr, False, False, True);
715 4 : tempInImageMask.copyData(*outmaskim_ptr);
716 : }
717 : else {
718 4 : tempInImageMask.copyData(inImageMask);
719 : }
720 :
721 : try {
722 8 : imregrid.regrid(chanMask, Interpolate2D::LINEAR, dirAxes, tempInImageMask);
723 0 : } catch (AipsError& x) {
724 0 : os<<LogIO::WARN<<"Regridding of the input mask image failed: "<<x.getMesg()<<LogIO::POST;
725 : }
726 : // extract input mask (first stokes plane)
727 16 : Array<Float> inMaskData;
728 16 : IPosition end2(4,outShape(outDirAxes(0)), outShape(outDirAxes(1)), 1, 1);
729 8 : if ( inNpol==outNpol ) {
730 4 : end2(outStokesAxis) = inNpol;
731 : }
732 8 : chanMask.doGetSlice(inMaskData, Slicer(start,end2));
733 16 : IPosition stride(4,1,1,1,1);
734 : //
735 : // continuum output mask case
736 8 : if (outNchan==1) {
737 : //No copying over channels, just do copying over all Stokes if input mask is a single Stokes
738 2 : if (inNpol == 1 && outNpol > 1) {
739 10 : for (Int ipol = 0; ipol < outNpol; ipol++) {
740 8 : start(outStokesAxis) = ipol;
741 8 : os<<"Copying input mask to Stokes plane="<<ipol<<LogIO::POST;
742 8 : outImageMask.putSlice(inMaskData,start,stride);
743 : }
744 : }
745 : }
746 : else { // for cube
747 24 : for (Int ich = 0; ich < outNchan; ich++) {
748 18 : start(outSpecAxis) = ich;
749 36 : IPosition inStart(4,0,0,0,0);
750 18 : if (inNpol == 1 && outNpol > 1) {
751 : // extend to other Stokes
752 30 : for (Int ipol = 0; ipol < outNpol; ipol++) {
753 24 : os<<"Copying input mask to Stokes plane="<<ipol<<LogIO::POST;
754 24 : start(outStokesAxis) = ipol;
755 24 : outImageMask.putSlice(inMaskData,start,stride);
756 6 : }
757 : }
758 : else {
759 : // copy Stokes plane as is (but expand it to all channels)
760 12 : stride(outStokesAxis) = stokesInc;
761 12 : if (inNpol == outNpol) {
762 42 : for (Int ipol = 0; ipol < outNpol; ipol++) {
763 : // need to slice mask from each stokes plane
764 30 : inMaskData.resize();
765 30 : inStart(outStokesAxis) = ipol;
766 30 : start(outStokesAxis) = ipol;
767 30 : end2(outStokesAxis) = 1;
768 30 : stride(outStokesAxis) = 1;
769 30 : chanMask.doGetSlice(inMaskData, Slicer(inStart,end2));
770 30 : outImageMask.putSlice(inMaskData,start,stride);
771 : }
772 :
773 : }
774 : else {
775 0 : outImageMask.putSlice(inMaskData,start,stride);
776 : }
777 : }
778 : } //for loop
779 : } // else
780 : }
781 : // Stokes I (a single Stokes plane mask with cube)
782 2 : else if (inNpol == 1) {
783 2 : if (inNpol != 1 ) {
784 0 : if (inShape(inStokesAxis)==outShape(outStokesAxis)) {
785 0 : stokesInc = inShape(inStokesAxis);
786 : }
787 : else {
788 0 : throw(AipsError("Cannot extend the input mask of "+String::toString(inNpol)+
789 0 : " Stokes planes to the output mask of "+String::toString(outNpol)+
790 0 : " Stokes planes. Please modify the input mask to make it a Stokes I mask or a mask with the same number of Stokes planes as the output.") );
791 : }
792 : }
793 2 : length(outStokesAxis) = stokesInc;
794 2 : length(outSpecAxis) = outNchan;
795 4 : Slicer slStokes(start, length);
796 : // make a subImage for regridding output (all channels)
797 4 : SubImage<Float> stokesMask(outImageMask, slStokes, true);
798 4 : ImageRegrid<Float> imregrid2;
799 6 : TempImage<Float> tempInStokesImageMask(stokesMask.shape(), stokesMask.coordinates(),memoryToUse());
800 2 : std::unique_ptr<ImageInterface<Float> > outstokesmaskim_ptr;
801 4 : Vector<Stokes::StokesTypes> outWhichPols;
802 2 : if ( stokesMask.shape().nelements() > inImageMask.shape().nelements() ) {
803 1 : casacore::ImageUtilities::addDegenerateAxes(os, outstokesmaskim_ptr, inImageMask,"",False, addSpecAxis, "I", False, False, True);
804 2 : Vector<Int> outWorldOrder(4);
805 2 : Vector<Int> outPixelOrder(4);
806 1 : outWorldOrder(0)=0;
807 1 : outWorldOrder(1)=1;
808 1 : outWorldOrder(2)=3;
809 1 : outWorldOrder(3)=2;
810 1 : outPixelOrder(0)=0;
811 1 : outPixelOrder(1)=1;
812 1 : outPixelOrder(2)=3;
813 1 : outPixelOrder(3)=2;
814 2 : CoordinateSystem modcsys=tempInStokesImageMask.coordinates();
815 2 : IPosition inMaskShape = inImageMask.shape();
816 2 : Array<Float> inData = outstokesmaskim_ptr->get();
817 2 : IPosition newAxisOrder(4,0,1,3,2);
818 2 : Array<Float> reorderedData = reorderArray(inData, newAxisOrder);
819 2 : IPosition newShape=reorderedData.shape();
820 : //os<< "reoderedData shape="<<reorderedData.shape()<<LogIO::POST;
821 2 : TempImage<Float> modTempInStokesMask(TiledShape(newShape), modcsys);
822 1 : modTempInStokesMask.put(reorderedData);
823 :
824 1 : if (compareSpectralCoordinate(inImageMask,tempInStokesImageMask) ) {
825 1 : tempInStokesImageMask.copyData(modTempInStokesMask);
826 : }
827 : }
828 : else {
829 1 : if (compareSpectralCoordinate(inImageMask,tempInStokesImageMask) ) {
830 1 : tempInStokesImageMask.copyData(inImageMask);
831 : }
832 :
833 : }
834 : try {
835 2 : imregrid2.regrid(stokesMask, Interpolate2D::LINEAR, dirAxes, tempInStokesImageMask);
836 0 : } catch (AipsError& x) {
837 0 : os<<LogIO::WARN<<"Regridding of the input mask image failed: "<<x.getMesg()<<LogIO::POST;
838 : }
839 2 : os <<"Slicing data..."<<LogIO::POST;
840 4 : Array<Float> inMaskData2;
841 4 : IPosition end3(4,outShape(outDirAxes(0)), outShape(outDirAxes(1)), 1, 1);
842 2 : end3(outStokesAxis) = inNpol;
843 2 : end3(outSpecAxis) = inNchan;
844 2 : stokesMask.doGetSlice(inMaskData2,slStokes);
845 4 : IPosition stride(4,1,1,1,1);
846 4 : IPosition inStart2(4,0,0,0,0);
847 10 : for (Int ipol = 0; ipol < outNpol; ipol++) {
848 : // need to slice mask from each stokes plane
849 8 : inMaskData2.resize();
850 8 : inStart2(outStokesAxis) = 0;
851 8 : start(outStokesAxis) = ipol;
852 8 : end3(outStokesAxis) = 1;
853 8 : stride(outStokesAxis) = 1;
854 8 : stride(outSpecAxis) = 1; // assume here inNchan == outNchan
855 8 : stokesMask.doGetSlice(inMaskData2, Slicer(inStart2,end3));
856 8 : outImageMask.putSlice(inMaskData2,start);
857 : }
858 :
859 : }
860 : else {
861 0 : throw(AipsError("Input mask,"+inImageMask.name()+" does not conform with the number of channels in output mask"));
862 : }
863 10 : }
864 :
865 2 : Bool SDMaskHandler::compareSpectralCoordinate(const ImageInterface<Float>& inImage, const ImageInterface<Float>& outImage)
866 : {
867 6 : LogIO os( LogOrigin("SDMaskHandler", "checkSpectralCoord",WHERE) );
868 :
869 4 : SpectralCoordinate outSpecCoord = outImage.coordinates().spectralCoordinate();
870 4 : IPosition inshape = inImage.shape();
871 4 : IPosition outshape = outImage.shape();
872 4 : CoordinateSystem incys = inImage.coordinates();
873 2 : CoordinateSystem outcsys = outImage.coordinates();
874 2 : Int inSpecAxis = CoordinateUtil::findSpectralAxis(incys);
875 2 : Int outSpecAxis = CoordinateUtil::findSpectralAxis(outcsys);
876 2 : Bool nchanMatch(true);
877 2 : if (inSpecAxis != -1 and outSpecAxis != -1 )
878 2 : nchanMatch = inshape(inSpecAxis) == outshape(outSpecAxis)? true: false;
879 2 : if (!nchanMatch) {
880 0 : if (!outSpecCoord.near(inImage.coordinates().spectralCoordinate())) {
881 0 : throw(AipsError("Cannot extend the input mask. Spectral coordiante and the number of channels of the input mask does not match with those of the output mask. Use a single channel mask or a mask that matches the spectral coordiante of the output. "));
882 : }
883 : else {
884 0 : throw(AipsError("Cannot extend the input mask. The number of the channels in Input mask,"+inImage.name()+"does not match with that of the output mask. Use a single channel mask or a mask that matches the spectral coordiante of the output. "));
885 : }
886 : return false;
887 : }
888 : else {
889 2 : if (!outSpecCoord.near(inImage.coordinates().spectralCoordinate())) {
890 0 : throw(AipsError("Cannot extend the input mask. Spectral coordiante of Input mask,"+inImage.name()+"does not match with that of the output mask. Use a single channel mask or a mask that matches the spectral coordiante of the output. "));
891 : return false;
892 : }
893 : }
894 2 : return true;
895 : }
896 :
897 :
898 : // was Imager::clone()...
899 : //static Bool cloneImShape(const ImageInterface<Float>& inImage, ImageInterface<Float>& outImage)
900 0 : Bool SDMaskHandler::cloneImShape(const ImageInterface<Float>& inImage, const String& outImageName)
901 : {
902 0 : LogIO os( LogOrigin("SDMaskHandler", "cloneImShape",WHERE) );
903 :
904 : try {
905 0 : PagedImage<Float> newImage(TiledShape(inImage.shape(),
906 0 : inImage.niceCursorShape()), inImage.coordinates(),
907 : // outImage.name());
908 0 : outImageName);
909 0 : newImage.set(0.0);
910 0 : newImage.table().flush(true, true);
911 0 : } catch (AipsError& x) {
912 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
913 0 : return false;
914 : }
915 0 : return true;
916 : }
917 :
918 :
919 0 : Int SDMaskHandler::makeInteractiveMask(std::shared_ptr<SIImageStore>& imstore,
920 : Int& niter, Int& cycleniter,
921 : String& threshold, String& cyclethreshold)
922 : {
923 : Int ret;
924 : // Int niter=1000, ncycles=100;
925 : // String thresh="0.001mJy";
926 0 : String imageName = imstore->getName()+".residual"+(imstore->getNTaylorTerms()>1?".tt0":"");
927 0 : String maskName = imstore->getName() + ".mask";
928 0 : imstore->mask()->unlock();
929 0 : cout << "Before interaction : niter : " << niter << " cycleniter : " << cycleniter << " thresh : " << threshold << " cyclethresh : " << cyclethreshold << endl;
930 : // ret = interactiveMasker_p->interactivemask(imageName, maskName,
931 : // niter, ncycles, threshold);
932 0 : return ret;
933 : }
934 :
935 0 : void SDMaskHandler::makeAutoMask(std::shared_ptr<SIImageStore> imstore)
936 : {
937 0 : LogIO os( LogOrigin("SDMaskHandler","makeAutoMask",WHERE) );
938 :
939 0 : Array<Float> localres;
940 : // Modification to be able to work with a cube (TT 2014-12-09)
941 : //imstore->residual()->get( localres , true );
942 0 : imstore->residual()->get( localres );
943 :
944 0 : Array<Float> localmask;
945 : //imstore->mask()->get( localmask , true );
946 0 : imstore->mask()->get( localmask );
947 :
948 0 : Int specAxis = CoordinateUtil::findSpectralAxis(imstore->mask()->coordinates());
949 0 : IPosition maskShape = localmask.shape();
950 0 : Int ndim = maskShape.nelements();
951 0 : IPosition pos(ndim,0);
952 0 : IPosition blc(ndim,0);
953 0 : IPosition trc(ndim,0);
954 0 : trc[0] = maskShape[0]-1;
955 0 : trc[1] = maskShape[1]-1;
956 : // added per channel mask setting
957 0 : for (pos[specAxis] = 0; pos[specAxis]<localmask.shape()[specAxis]; pos[specAxis]++)
958 : {
959 0 : IPosition posMaxAbs( localmask.shape().nelements(), 0);
960 0 : blc[specAxis]=pos[specAxis];
961 0 : trc[specAxis]=pos[specAxis];
962 0 : Float maxAbs=0.0;
963 : Float minVal;
964 0 : IPosition posmin(localmask.shape().nelements(), 0);
965 : //minMax(minVal, maxAbs, posmin, posMaxAbs, localres);
966 0 : minMax(minVal, maxAbs, posmin, posMaxAbs, localres(blc,trc));
967 :
968 : // cout << "Max position : " << posMaxAbs << endl;
969 :
970 0 : Int dist=5;
971 :
972 : //IPosition pos(2,0,0); // Deal with the input shapes properly......
973 0 : for (pos[0]=posMaxAbs[0]-dist; pos[0]<posMaxAbs[0]+dist; pos[0]++)
974 : {
975 0 : for (pos[1]=posMaxAbs[1]-dist; pos[1]<posMaxAbs[1]+dist; pos[1]++)
976 : {
977 0 : if( pos[0]>0 && pos[0]<localmask.shape()[0] && pos[1]>0 && pos[1]<localmask.shape()[1] )
978 : {
979 0 : localmask( pos ) = 1.0;
980 : }
981 : }
982 : }
983 : } // over channels
984 : //cout << "Sum of mask : " << sum(localmask) << endl;
985 0 : Float summask = sum(localmask);
986 0 : if( summask==0.0 ) { localmask=1.0; summask = sum(localmask); }
987 0 : os << LogIO::NORMAL1 << "Make Autobox mask with " << summask << " available pixels " << LogIO::POST;
988 :
989 0 : imstore->mask()->put( localmask );
990 :
991 : // imstore->mask()->get( localmask , true );
992 : // cout << "Sum of imstore mask : " << sum( localmask ) << endl;
993 :
994 0 : }
995 :
996 52 : void SDMaskHandler::autoMask(std::shared_ptr<SIImageStore> imstore,
997 : ImageInterface<Float>& posmask,
998 : const Int iterdone,
999 : Vector<Bool>& chanflag,
1000 : Record& robuststatsrec,
1001 : const String& alg,
1002 : const String& threshold,
1003 : const Float& fracofpeak,
1004 : const String& resolution,
1005 : const Float& resbybeam,
1006 : const Int nmask,
1007 : const Bool autoadjust,
1008 : // new params for the multithreshold alg.
1009 : const Float& sidelobethreshold,
1010 : const Float& noisethreshold,
1011 : const Float& lownoisethreshold,
1012 : const Float& negativethreshold,
1013 : const Float& cutthreshold,
1014 : const Float& smoothfactor,
1015 : const Float& minbeamfrac,
1016 : const Int growiterations,
1017 : const Bool dogrowprune,
1018 : const Float& minpercentchange,
1019 : const Bool verbose,
1020 : const Bool fastnoise,
1021 : const Bool isthresholdreached,
1022 : Float pblimit)
1023 : {
1024 104 : LogIO os( LogOrigin("SDMaskHandler","autoMask",WHERE) );
1025 :
1026 : //currently supported alg:
1027 : // onebox: a box around a max (box size +/-5pix around max position)
1028 : // thresh: threshold based auto masking (uses threshold or fracofpeak, and resolution)
1029 : //
1030 : // create a working copy of residual image (including pixel masks)
1031 52 : os << LogIO::NORMAL2 <<"algorithm:"<<alg<<LogIO::POST;
1032 52 : TempImage<Float>* tempres = new TempImage<Float>(imstore->residual()->shape(), imstore->residual()->coordinates(), memoryToUse());
1033 : //Array<Float> maskdata;
1034 : //Array<Float> psfdata;
1035 : {
1036 52 : Array<Float> resdata;
1037 52 : imstore->residual()->get(resdata);
1038 52 : tempres->put(resdata);
1039 52 : tempres->setImageInfo(imstore->residual()->imageInfo());
1040 52 : tempres->attachMask(ArrayLattice<Bool> (imstore->residual()->getMask()));
1041 : }
1042 :
1043 52 : TempImage<Float>* tempmask = new TempImage<Float>(imstore->mask()->shape(), imstore->mask()->coordinates(), memoryToUse());
1044 : // get current mask and apply pbmask
1045 : {
1046 104 : Array<Float> maskdata;
1047 52 : imstore->mask()->get(maskdata);
1048 104 : String maskname = imstore->getName()+".mask";
1049 52 : tempmask->put(maskdata);
1050 : //For debug
1051 : //PagedImage<Float> tempmaskInit(tempmask->shape(), tempmask->coordinates(), "tempMaskInit"+String::toString(iterdone));
1052 : //tempmaskInit.copyData(*tempmask);
1053 :
1054 : //
1055 :
1056 52 : if (pblimit>0.0 && imstore->hasPB()) {
1057 : //Determine if there is mask and set pixel mask to True for the mask to check by ntrue
1058 0 : LatticeExpr<Bool> pixmask( iif(*tempmask > 0.0, True, False));
1059 0 : TempImage<Float>* dummy = new TempImage<Float>(tempres->shape(), tempres->coordinates(), memoryToUse());
1060 0 : dummy->attachMask(pixmask);
1061 0 : LatticeExpr<Float> themask;
1062 0 : if (!ntrue(dummy->getMask())) { // initial zero mask
1063 : //os<<"INITIAL Zero Mask ...."<<LogIO::POST;
1064 : //themask = LatticeExpr<Float>( iif( (*(imstore->pb())) > pblimit , 1.0 , 0.0 ));
1065 0 : themask = LatticeExpr<Float>( *tempmask);
1066 : }
1067 : else {
1068 : //os<<"INITIAL NON-Zero Mask ...."<<LogIO::POST;
1069 0 : themask = LatticeExpr<Float>( iif( (*(imstore->pb())) > pblimit, *(imstore->mask()), 0.0));
1070 : }
1071 : // attache pixmask to temp res image to be used in stats etc
1072 0 : tempres->attachMask(LatticeExpr<Bool> ( iif(*(imstore->pb()) > pblimit, True, False)));
1073 0 : imstore->mask()->copyData( themask );
1074 0 : imstore->mask()->get(maskdata);
1075 0 : tempmask->put(maskdata);
1076 0 : delete dummy;
1077 : }
1078 : }
1079 : //for debug
1080 : //String tempresname="initialRes_"+String::toString(iterdone)+".im";
1081 : //PagedImage<Float> initialRes(tempres->shape(), tempres->coordinates(), tempresname);
1082 : //initialRes.copyData(*tempres);
1083 :
1084 : // Not use this way for now. Got an issue on removing pixel mask from *.mask image
1085 : // retrieve pixelmask (i.e. pb mask)
1086 : //LatticeExpr<Bool> pixmasyyk;
1087 : //if (imstore->mask()->hasPixelMask()) {
1088 : // pixmask = LatticeExpr<Bool> (imstore->mask()->pixelMask());
1089 : //
1090 : // create pixel mask (set to True for the previous selected region(s) to exclude the region from the stats/masking )
1091 : //LatticeExpr<Bool> prevmask( iif(*tempmask > 0.0 || pixmask, True, False) );
1092 : // TempImage<Float>* dummy = new TempImage<Float>(tempres->shape(), tempres->coordinates());
1093 : // dummy->attachMask(pixmask);
1094 : //if (ntrue(dummy->getMask())) tempres->attachMask(pixmask);
1095 : // if (ntrue(dummy->getMask())) {
1096 : //tempmask->removeMask();
1097 : // }
1098 : // else {
1099 : // os<<LogIO::DEBUG1<<"No pixel mask"<<LogIO::POST;
1100 : // }
1101 : // delete dummy; dummy=0;
1102 : //}
1103 : //
1104 : //input
1105 52 : Quantity qthresh(0,"");
1106 52 : Quantity qreso(0,"");
1107 52 : Quantity::read(qreso,resolution);
1108 52 : Float sigma = 0.0;
1109 : // if fracofpeak (fraction of a peak) is specified, use it to set a threshold
1110 52 : if ( fracofpeak != 0.0 ) {
1111 0 : if (fracofpeak > 1.0 )
1112 0 : throw(AipsError("Fracofpeak must be < 1.0"));
1113 0 : sigma = 0.0;
1114 : }
1115 52 : else if(Quantity::read(qthresh,threshold) ) {
1116 : // evaluate threshold input
1117 : //cerr<<"qthresh="<<qthresh.get().getValue()<<" unit="<<qthresh.getUnit()<<endl;
1118 52 : if (qthresh.getUnit()!="") {
1119 : // use qthresh and set sigma =0.0 to ignore
1120 0 : sigma = 0.0;
1121 : }
1122 : else {
1123 52 : sigma = String::toFloat(threshold);
1124 52 : if (sigma==0.0) {
1125 : // default case: threshold, fracofpeak unset => use default (3sigma)
1126 52 : sigma = 3.0;
1127 : }
1128 : }
1129 : }
1130 : else {
1131 0 : if (!sigma) {
1132 0 : throw(AipsError("Unrecognized automask threshold="+threshold));
1133 : }
1134 : }
1135 :
1136 : //TempImage<Float>* resforstats = new TempImage<Float>(imstore->residual()->shape(), imstore->residual()->coordinates());
1137 : //Array<Float> resdata2;
1138 : //imstore->residual()->get(resdata2);
1139 : //resforstats->put(resdata2);
1140 : //resforstats->setImageInfo(imstore->residual()->imageInfo());
1141 : //LatticeExpr<Bool> prevmask( iif(*tempmask > 0.0 , True, False) );
1142 : //resforstats->attachMask(prevmask);
1143 : //std::shared_ptr<casacore::ImageInterface<float> > tempres_ptr(resforstats);
1144 : /**
1145 : std::shared_ptr<casacore::ImageInterface<float> > tempres_ptr(tempres);
1146 : //cerr<<" tempres->hasPixelMask? "<<tempres->hasPixelMask()<<endl;
1147 : ImageStatsCalculator imcalc( tempres_ptr, 0, "", False);
1148 : Vector<Int> axes(2);
1149 : axes[0] = 0;
1150 : axes[1] = 1;
1151 : imcalc.setAxes(axes);
1152 : // for now just for new autobox alg.
1153 : if (alg.contains("newauto") ) {
1154 : imcalc.setRobust(true);
1155 : }
1156 : Record thestats = imcalc.statistics();
1157 : ***/
1158 :
1159 52 : Record* region_ptr=0;
1160 52 : String LELmask("");
1161 : // note: tempres is res image with pblimit applied
1162 52 : Bool robust(false); // robust is false by default (turn it on for 'multithresh')
1163 52 : Bool doAnnulus(false); // turn off stats on an annulus
1164 52 : if (alg.contains("multithresh") ) {
1165 52 : robust=true;
1166 : // define an annulus
1167 52 : Float outerRadius=0.0;
1168 52 : Float innerRadius=1.0;
1169 52 : if (imstore->hasPB()) {
1170 43 : LatticeExpr<Bool> blat;
1171 : //use pb based annulus for statistics
1172 43 : if (doAnnulus) {
1173 0 : outerRadius = 0.2;
1174 0 : innerRadius = 0.3;
1175 0 : blat=LatticeExpr<Bool> (iif(( *imstore->pb() < innerRadius && *imstore->pb() > outerRadius) , True, False) );
1176 : }
1177 : else {
1178 43 : blat=LatticeExpr<Bool> ((*imstore->pb()).pixelMask());
1179 : }
1180 43 : TempImage<Float>* testres = new TempImage<Float>(imstore->residual()->shape(), imstore->residual()->coordinates(), memoryToUse());
1181 43 : testres->set(1.0);
1182 43 : testres->attachMask(blat);
1183 :
1184 43 : if (ntrue(testres->getMask()) > 0 ) {
1185 86 : String pbname = imstore->getName()+".pb";
1186 43 : if (doAnnulus) {
1187 0 : LELmask=pbname+"<"+String::toString(innerRadius)+" && "+pbname+">"+String::toString(outerRadius);
1188 : }
1189 : }
1190 43 : delete testres; testres=0;
1191 : }
1192 : }
1193 :
1194 : //use new noise calc.
1195 : //Bool useoldstats(False);
1196 :
1197 : // at this point tempres has pbmask applied
1198 52 : LatticeExpr<Bool> pbmask(tempres->pixelMask());
1199 :
1200 52 : Record thestats = calcImageStatistics(*tempres, LELmask, region_ptr, robust);
1201 52 : Array<Double> maxs, mins, rmss, mads;
1202 52 : thestats.get(RecordFieldId("max"), maxs);
1203 52 : thestats.get(RecordFieldId("rms"), rmss);
1204 :
1205 52 : Record thenewstats;
1206 52 : if (!(iterdone==0 && robuststatsrec.nfields()) ) { // this is an indirect way to check if initial stats by nsigma threshold is already run.
1207 48 : if (fastnoise) { // use a faster but less robust noise determination
1208 44 : thenewstats = thestats;
1209 44 : os<< LogIO::DEBUG1 << " *** Using classic image statistics ***"<<LogIO::POST;
1210 44 : os<< LogIO::DEBUG1 << "All rms's on the input image -- rms.nelements()="<<rmss.nelements()<<" rms="<<rmss<<LogIO::POST;
1211 : //os<< LogIO::DEBUG1 << "All max's on the input image -- max.nelements()="<<maxs.nelements()<<" max="<<maxs<<LogIO::POST;
1212 44 : if (alg.contains("multithresh")) {
1213 44 : thestats.get(RecordFieldId("medabsdevmed"), mads);
1214 44 : os<< LogIO::DEBUG1 << "All MAD's on the input image -- mad.nelements()="<<mads.nelements()<<" mad="<<mads<<LogIO::POST;
1215 : }
1216 : }
1217 : else {
1218 : // Revised version of calcRobustImageStatistics (previous one- rename to calcRobustImageStatisticsOld)
1219 : //Record thenewstats = calcRobustImageStatistics(*tempres, *tempmask, pbmask, LELmask, region_ptr, robust, chanflag);
1220 : try{
1221 4 : thenewstats = calcRobustImageStatistics(*tempres, *tempmask, pbmask, LELmask, region_ptr, robust, chanflag);
1222 : }
1223 0 : catch( AipsError &x )
1224 : {
1225 : //now that there are part images that are masked...should just not proceed rather than throw an exception
1226 0 : if(x.getMesg().contains("zero element data"))
1227 0 : return;
1228 : else
1229 0 : throw(x);
1230 : }
1231 8 : Array<Double> newmaxs, newmins, newrmss, newmads;
1232 4 : thenewstats.get(RecordFieldId("max"), newmaxs);
1233 4 : thenewstats.get(RecordFieldId("rms"), newrmss);
1234 4 : os<< LogIO::DEBUG1 << "*** Using the new image statistics *** "<<LogIO::POST;
1235 4 : os<< LogIO::DEBUG1 << "All NEW rms's on the input image -- rms.nelements()="<<newrmss.nelements()<<" rms="<<newrmss<<LogIO::POST;
1236 : //os<< LogIO::DEBUG1 << "All NEW max's on the input image -- max.nelements()="<<newmaxs.nelements()<<" max="<<newmaxs<<LogIO::POST;
1237 4 : if (alg.contains("multithresh")) {
1238 4 : thenewstats.get(RecordFieldId("medabsdevmed"), newmads);
1239 4 : os<< LogIO::DEBUG1 << "All NEW MAD's on the input image -- mad.nelements()="<<newmads.nelements()<<" mad="<<newmads<<LogIO::POST;
1240 : }
1241 : }
1242 48 : os<<LogIO::NORMAL3<<" set the stats to what is calculated here" <<LogIO::POST;
1243 48 : robuststatsrec=thenewstats;
1244 : }
1245 : else {
1246 4 : thenewstats=robuststatsrec; // use existing
1247 : }
1248 : //os<<" current robuststatsrec nfields="<<robuststatsrec.nfields()<<LogIO::POST;
1249 :
1250 :
1251 52 : os<<LogIO::NORMAL <<"SidelobeLevel = "<<imstore->getPSFSidelobeLevel()<<LogIO::POST;
1252 52 : itsSidelobeLevel = imstore->getPSFSidelobeLevel();
1253 : //os<< "mask algortihm ="<<alg<< LogIO::POST;
1254 52 : if (alg==String("") || alg==String("onebox")) {
1255 : //cerr<<" calling makeAutoMask(), simple 1 cleanbox around the max"<<endl;
1256 0 : makeAutoMask(imstore);
1257 : }
1258 52 : else if (alg==String("thresh")) {
1259 0 : autoMaskByThreshold(*tempmask, *tempres, *imstore->psf(), qreso, resbybeam, qthresh, fracofpeak,
1260 : thestats, sigma, nmask, autoadjust);
1261 : }
1262 52 : else if (alg==String("thresh2")) {
1263 0 : autoMaskByThreshold2(*tempmask, *tempres, *imstore->psf(), qreso, resbybeam, qthresh, fracofpeak, thestats, sigma, nmask);
1264 : }
1265 52 : else if (alg==String("multithresh")) {
1266 52 : autoMaskByMultiThreshold(*tempmask, posmask, *tempres, *imstore->psf(), thestats, thenewstats, iterdone, chanflag, minpercentchange, itsSidelobeLevel, sidelobethreshold, noisethreshold, lownoisethreshold, negativethreshold, cutthreshold, smoothfactor, minbeamfrac, growiterations, dogrowprune, verbose, isthresholdreached);
1267 : }
1268 :
1269 : // this did not work (it won't physically remove the mask from the image
1270 : /***
1271 : if (imstore->mask()->hasPixelMask()) {
1272 : imstore.get()->mask()->removeRegion(fname, RegionHandler::Any, False);
1273 : cerr<<"imstore->mask()->name()="<<imstore->mask()->name()<<endl;
1274 : cerr<<" mask: "<<fname<<" exist on disk? ="<<File(fname).isDirectory()<<endl;
1275 : }
1276 : ***/
1277 : {
1278 52 : Array<Float> updatedMaskData;
1279 52 : tempmask->get(updatedMaskData);
1280 52 : imstore->mask()->put(updatedMaskData);
1281 : }
1282 : //tempmask->get(maskdata);
1283 : //imstore->mask()->put(maskdata);
1284 52 : delete tempmask; tempmask=0;
1285 : //delete temppsf; temppsf=0;
1286 52 : delete tempres; tempres=0;
1287 : }
1288 :
1289 743 : Record SDMaskHandler::calcImageStatistics(ImageInterface<Float>& res, String& LELmask, Record* regionPtr, const Bool robust )
1290 : {
1291 2229 : LogIO os( LogOrigin("SDMaskHandler","calcImageStatistics",WHERE) );
1292 743 : TempImage<Float>* tempres = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
1293 1486 : Array<Float> resdata;
1294 : //
1295 :
1296 743 : res.get(resdata);
1297 743 : tempres->put(resdata);
1298 : // if input image (res) has a pixel mask, make sure to honor it so the region is exclude from statistics
1299 743 : if (res.hasPixelMask()) {
1300 163 : tempres->attachMask(res.pixelMask());
1301 : }
1302 :
1303 743 : tempres->setImageInfo(res.imageInfo());
1304 1486 : std::shared_ptr<casacore::ImageInterface<Float> > tempres_ptr(tempres);
1305 :
1306 : // 2nd arg is regionRecord, 3rd is LELmask expression and those will be AND
1307 : // to define a region to be get statistics
1308 : //ImageStatsCalculator imcalc( tempres_ptr, 0, "", False);
1309 : //String lelstring = pbname+">0.92 && "+pbname+"<0.98";
1310 : //cerr<<"lelstring = "<<lelstring<<endl;
1311 : //cerr<<"LELMask="<<LELmask<<endl;
1312 1486 : ImageStatsCalculator<Float> imcalc( tempres_ptr, regionPtr, LELmask, False);
1313 1486 : Vector<Int> axes(2);
1314 743 : axes[0] = 0;
1315 743 : axes[1] = 1;
1316 743 : imcalc.setAxes(axes);
1317 743 : imcalc.setRobust(robust);
1318 743 : Record thestats = imcalc.statistics();
1319 :
1320 : //cerr<<"thestats="<<thestats<<endl;
1321 : //Array<Double> max, min, rms, mad;
1322 : //thestats.get(RecordFieldId("max"), max);
1323 1486 : return thestats;
1324 : }
1325 :
1326 : // robust image statistics for better noise estimation
1327 0 : Record SDMaskHandler::calcRobustImageStatisticsOld(ImageInterface<Float>& res, ImageInterface<Float>& prevmask , LatticeExpr<Bool>& pbmask, String& LELmask, Record* regionPtr, const Bool robust, Vector<Bool>& chanflag )
1328 : {
1329 0 : LogIO os( LogOrigin("SDMaskHandler","calcRobustImageStatisticsOld",WHERE) );
1330 :
1331 : //Bool debugStats(true);
1332 :
1333 0 : Array<Float> wholemaskdata;
1334 0 : IPosition maskshp = prevmask.shape();
1335 0 : IPosition maskstart(4,0);
1336 0 : IPosition masklength(4,maskshp(0),maskshp(1), 1, 1);
1337 0 : Slicer masksl(maskstart, masklength);
1338 0 : prevmask.doGetSlice(wholemaskdata,masksl); // single plane
1339 0 : Float npixwholemask = sum(wholemaskdata);
1340 0 : Bool fullmask(False);
1341 0 : if (npixwholemask == (Float) maskshp(0) * (Float) maskshp(1) ) {
1342 0 : fullmask=True;
1343 0 : os<<LogIO::DEBUG1 <<"Appears to be fully masked! npix for the whole mask ="<<npixwholemask<<LogIO::POST;
1344 : }
1345 : //TempImage<Float>* tempres = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
1346 : //Array<Float> resdata;
1347 :
1348 : //res.get(resdata);
1349 : //tempres->put(resdata);
1350 : // if input image (res) has a pixel mask, make sure to honor it so the region is exclude from statistics
1351 : //if (res.hasPixelMask()) {
1352 : // tempres->attachMask(res.pixelMask());
1353 : //}
1354 0 : if (nfalse(pbmask.getMask())) {
1355 0 : os<<LogIO::DEBUG1<<"has pbmask"<<LogIO::POST;
1356 : }
1357 :
1358 : //check Stokes and spectral axes
1359 : //IPosition shp = tempres->shape();
1360 0 : IPosition shp = res.shape();
1361 : //Int specaxis = CoordinateUtil::findSpectralAxis(tempres->coordinates());
1362 0 : Int specaxis = CoordinateUtil::findSpectralAxis(res.coordinates());
1363 0 : uInt nchan = shp(specaxis);
1364 0 : Vector<Stokes::StokesTypes> whichPols;
1365 : //Int stokesaxis = CoordinateUtil::findStokesAxis(whichPols, tempres->coordinates());
1366 0 : Int stokesaxis = CoordinateUtil::findStokesAxis(whichPols, res.coordinates());
1367 0 : uInt nStokes= shp(stokesaxis);
1368 : uInt iaxis2;
1369 : uInt iaxis3;
1370 :
1371 : //TempImage<Bool> pbmaskim(shp, tempres->coordinates(), memoryToUse());
1372 0 : TempImage<Bool> pbmaskim(shp, res.coordinates(), memoryToUse());
1373 0 : pbmaskim.copyData(pbmask);
1374 :
1375 : //check dimensions of tempres, prevmask, pbmask
1376 0 : if (prevmask.shape()!=shp) {
1377 0 : throw(AipsError("Mismatch in shapes of the previous mask and residual image"));
1378 : }
1379 0 : else if (pbmask.shape()!=shp) {
1380 0 : throw(AipsError("Mismatch in shapes of pbmask and residual image"));
1381 : }
1382 :
1383 : //for stats storage
1384 0 : IPosition statshape(2, shp(2), shp(3));
1385 0 : Array<Double> outMins(statshape);
1386 0 : Array<Double> outMaxs(statshape);
1387 0 : Array<Double> outRmss(statshape);
1388 0 : Array<Double> outMads(statshape);
1389 0 : Array<Double> outMdns(statshape);
1390 :
1391 0 : for (uInt istokes = 0; istokes < nStokes; istokes++) {
1392 0 : std::vector<double> mins, maxs, rmss, mads, mdns;
1393 0 : for (uInt ichan = 0; ichan < nchan; ichan++ ) {
1394 0 : if (stokesaxis==3) {
1395 0 : iaxis2 = ichan;
1396 0 : iaxis3 = istokes;
1397 : }
1398 : else {
1399 0 : iaxis2 = istokes;
1400 0 : iaxis3 = ichan;
1401 : }
1402 0 : if (chanflag.nelements()==0 || !chanflag(ichan)) { // get new stats
1403 : // check if mask is empty (evaulated as the whole input mask )
1404 0 : Array<Float> maskdata;
1405 : //IPosition maskshape = prevmask.shape();
1406 : //Int naxis = maskshape.size();
1407 : //IPosition blc(naxis,0);
1408 : //IPosition trc=maskshape-1;
1409 : //Slicer sl(blc,trc,Slicer::endIsLast);
1410 0 : IPosition start(4, 0, 0, iaxis2, iaxis3);
1411 0 : IPosition length(4, shp(0),shp(1), 1, 1);
1412 0 : Slicer sl(start, length);
1413 : //os<<"slicer for subimage start="<<start<<" length="<<length<<LogIO::POST;
1414 : // create subimage (slice) of tempres
1415 0 : AxesSpecifier aspec(False);
1416 0 : SubImage<Float>* subRes = new SubImage<Float>(res, sl, aspec, True);
1417 0 : TempImage<Float>* tempSubRes = new TempImage<Float>(subRes->shape(), subRes->coordinates(), memoryToUse());
1418 0 : tempSubRes->copyData(*subRes);
1419 0 : SubImage<Float>* subprevmask = new SubImage<Float>(prevmask, sl, aspec, True);
1420 0 : SubImage<Bool>* subpbmask = new SubImage<Bool>(pbmaskim, sl, aspec, True);
1421 :
1422 0 : prevmask.doGetSlice(maskdata,sl); // single plane
1423 0 : Float nmaskpix = sum(maskdata);
1424 0 : Array<Bool> booldata;
1425 0 : pbmask.doGetSlice(booldata,sl); // single plane pbmask
1426 : //os<<"valid pix ntrue(booldata)="<<ntrue(booldata)<<LogIO::POST;
1427 :
1428 0 : if (nmaskpix==0 ) {
1429 0 : if (ntrue(booldata)) {
1430 0 : os<<LogIO::DEBUG1<<"No existing mask but apply pbmask..."<<LogIO::POST;
1431 0 : tempSubRes->attachMask(*subpbmask);
1432 : //os<<"ntrue tempres pixelmask=="<<ntrue(tempres->getMask())<<LogIO::POST;
1433 : }
1434 : }
1435 0 : else if (!fullmask) {
1436 0 : os<<LogIO::DEBUG1<<"Do the image statitistics in a region outside the mask..."<<LogIO::POST;
1437 0 : LatticeExpr<Bool> outsideMaskReg;
1438 : //if (res.hasPixelMask()) {
1439 0 : if (nfalse(booldata)) {
1440 : //LatticeExpr<Bool> pbmask(iif(res.pixelMask());
1441 0 : outsideMaskReg = LatticeExpr<Bool> (iif(*subprevmask == 1.0 || !*subpbmask, False, True));
1442 : // for debug
1443 : //if (debugStats) {
1444 : // PagedImage<Float> pbmaskSave(res.shape(), res.coordinates(), "pbmasksaved.im");
1445 : // LatticeExpr<Float> temppbmasksave(iif(subpbmask, 1.0, 0.0));
1446 : // pbmaskSave.copyData(temppbmasksave);
1447 : //}
1448 : }
1449 : else {
1450 0 : outsideMaskReg = LatticeExpr<Bool> (iif(*subprevmask == 1.0, False, True));
1451 : //LatticeExpr<Bool> outsideMaskReg(iif((*subprevmask == 1.0 || !*subpbmask, False, True));
1452 : }
1453 : //tempres->attachMask(outsideMaskReg);
1454 0 : tempSubRes->attachMask(outsideMaskReg);
1455 : }
1456 :
1457 : //DEBUG
1458 : //if(debugStats) {
1459 : // PagedImage<Float> temptempIm(res.shape(), res.coordinates(), "temptempmask.im");
1460 : // temptempIm.copyData(prevmask);
1461 : // PagedImage<Float> temptempResIm(res.shape(), res.coordinates(), "temptempres.im");
1462 : // temptempResIm.copyData(*tempres);
1463 : // temptempResIm.makeMask("maskcopy",True, True);
1464 : // if (tempres->hasPixelMask()) {
1465 : // temptempResIm.pixelMask().put((tempres->pixelMask()).get());
1466 : // }
1467 : //} //for debug
1468 :
1469 0 : std::shared_ptr<casacore::ImageInterface<Float> > tempres_ptr(tempSubRes);
1470 :
1471 : // 2nd arg is regionRecord, 3rd is LELmask expression and those will be AND
1472 : // to define a region to be get statistics
1473 : //ImageStatsCalculator imcalc( tempres_ptr, 0, "", False);
1474 0 : ImageStatsCalculator<Float> imcalc( tempres_ptr, regionPtr, LELmask, True);
1475 0 : Vector<Int> axes(2);
1476 0 : axes[0] = 0;
1477 0 : axes[1] = 1;
1478 0 : imcalc.setAxes(axes);
1479 :
1480 : // for an empty (= no mask) mask, use Chauvenet algorithm with maxiter=5 and zscore=-1
1481 0 : if (nmaskpix==0.0 || fullmask ) {
1482 0 : os<<"Using Chauvenet algorithm for image statistics"<<LogIO::POST;
1483 0 : imcalc.configureChauvenet(Double(-1.0), Int(5));
1484 : }
1485 0 : imcalc.setRobust(robust);
1486 0 : Record thestats = imcalc.statistics();
1487 :
1488 0 : Array<Double> arrmins, arrmaxs, arrrmss, arrmads, arrmdns;
1489 :
1490 : // repack
1491 0 : thestats.get(RecordFieldId("min"), arrmins);
1492 0 : thestats.get(RecordFieldId("max"), arrmaxs);
1493 0 : thestats.get(RecordFieldId("rms"), arrrmss);
1494 : //robust = True only
1495 0 : if (robust) {
1496 0 : thestats.get(RecordFieldId("medabsdevmed"), arrmads);
1497 0 : thestats.get(RecordFieldId("median"), arrmdns);
1498 : }
1499 0 : if (arrmaxs.nelements()==0 ){
1500 0 : throw(AipsError("No image statisitics is returned. Possible the whole image is masked."));
1501 : }
1502 0 : IPosition zeroindx(arrmins.ndim(), 0);
1503 :
1504 0 : mins.push_back(arrmins(zeroindx));
1505 0 : maxs.push_back(arrmaxs(zeroindx));
1506 0 : rmss.push_back(arrrmss(zeroindx));
1507 0 : if (robust) {
1508 0 : mads.push_back(arrmads(zeroindx));
1509 0 : mdns.push_back(arrmdns(zeroindx));
1510 : }
1511 : //os<<"deleting tempSubRes"<<LogIO::POST;
1512 : //delete tempSubRes; tempSubRes=0;
1513 0 : delete subRes; subRes=0;
1514 0 : delete subprevmask; subprevmask=0;
1515 0 : delete subpbmask; subpbmask=0;
1516 :
1517 : }// if-ichanflag end
1518 : else {
1519 0 : mins.push_back(Double(0.0));
1520 0 : maxs.push_back(Double(0.0));
1521 0 : rmss.push_back(Double(0.0));
1522 0 : if (robust) {
1523 0 : mads.push_back(Double(0.0));
1524 0 : mdns.push_back(Double(0.0));
1525 : }
1526 : }
1527 : } // chan-for-loop end
1528 : //os<<" rms vector for stokes="<<istokes<<" : "<<rmss<<LogIO::POST;
1529 : //os<<"outMins.shape="<<outMins.shape()<<LogIO::POST;
1530 :
1531 0 : IPosition start(2, istokes, 0);
1532 0 : IPosition length(2, 1, nchan); // slicer per Stokes
1533 : //Slicer sl(blc, trc,Slicer::endIsLast );
1534 0 : Slicer slstokes(start, length);
1535 : //cerr<<"set outMin slstokes="<<slstokes<<" to the mins vector"<<endl;
1536 0 : Vector<Double> minvec(mins);
1537 0 : Vector<Double> maxvec(maxs);
1538 0 : Vector<Double> rmsvec(rmss);
1539 : //os<<"stl vector rmss"<<rmss<<LogIO::POST;
1540 : //os<<"rmsvec.shape="<<rmsvec.shape()<<" rmsvec="<<rmsvec<<LogIO::POST;
1541 0 : Matrix<Double> minmat(minvec);
1542 0 : Matrix<Double> maxmat(maxvec);
1543 : //Matrix<Double> rmsmat((Vector<Double>) rmss);
1544 0 : Matrix<Double> rmsmat(rmsvec);
1545 : //os<<"Intial shape of rmsmat="<<rmsmat.shape()<<LogIO::POST;
1546 : //os<<"Intial content of rmsmat="<<rmsmat<<LogIO::POST;
1547 : // copyValues do not work if there is a degerate axis in front
1548 0 : minmat.resize(IPosition(2, nchan, 1),True);
1549 0 : maxmat.resize(IPosition(2, nchan, 1),True);
1550 0 : rmsmat.resize(IPosition(2, nchan, 1),True);
1551 0 : Array<Double> minarr = transpose(minmat);
1552 0 : Array<Double> maxarr = transpose(maxmat);
1553 0 : Array<Double> rmsarr = transpose(rmsmat);
1554 0 : if (robust) {
1555 0 : Vector<Double> madvec(mads);
1556 0 : Vector<Double> mdnvec(mdns);
1557 0 : Matrix<Double> madmat(madvec);
1558 0 : Matrix<Double> mdnmat(mdnvec);
1559 0 : madmat.resize(IPosition(2, nchan, 1),True);
1560 0 : mdnmat.resize(IPosition(2, nchan, 1),True);
1561 0 : Array<Double> madarr = transpose(madmat);
1562 0 : Array<Double> mdnarr = transpose(mdnmat);
1563 0 : outMads(slstokes) = madarr;
1564 0 : outMdns(slstokes) = mdnarr;
1565 : //os<<"madarr="<<madarr<<LogIO::POST;
1566 : //os<<"mdnarr="<<mdnarr<<LogIO::POST;
1567 : }
1568 : //os<<"rmsarr="<<rmsarr<<LogIO::POST;
1569 0 : outMins(slstokes) = minarr;
1570 0 : outMaxs(slstokes) = maxarr;
1571 0 : outRmss(slstokes) = rmsarr;
1572 : //cerr<<"Done setting outMin sltokes="<<slstokes<<" to the mins vector"<<endl;
1573 : } // stokes-for-loop end
1574 :
1575 0 : Record theOutStatRec;
1576 0 : if (shp(2) == 1) { //Single Stokes plane
1577 : // remove degenerate axis
1578 0 : theOutStatRec.define("min", outMins.nonDegenerate(IPosition(1,1)));
1579 0 : theOutStatRec.define("max", outMaxs.nonDegenerate(IPosition(1,1)));
1580 0 : theOutStatRec.define("rms", outRmss.nonDegenerate(IPosition(1,1)));
1581 0 : if (robust) {
1582 0 : theOutStatRec.define("medabsdevmed", outMads.nonDegenerate(IPosition(1,1)));
1583 0 : theOutStatRec.define("median", outMdns.nonDegenerate(IPosition(1,1)));
1584 : }
1585 : }
1586 0 : else if(shp(3) == 1) { //Single channel plane
1587 0 : theOutStatRec.define("min", outMins.nonDegenerate(IPosition(1,0)));
1588 0 : theOutStatRec.define("max", outMaxs.nonDegenerate(IPosition(1,0)));
1589 0 : theOutStatRec.define("rms", outRmss.nonDegenerate(IPosition(1,0)));
1590 0 : if (robust) {
1591 0 : theOutStatRec.define("medabsdevmed", outMads.nonDegenerate(IPosition(1,0)));
1592 0 : theOutStatRec.define("median", outMdns.nonDegenerate(IPosition(1,0)));
1593 : }
1594 : }
1595 : else {
1596 0 : theOutStatRec.define("min", outMins);
1597 0 : theOutStatRec.define("max", outMaxs);
1598 0 : theOutStatRec.define("rms", outRmss);
1599 0 : if (robust) {
1600 0 : theOutStatRec.define("medabsdevmed", outMads);
1601 0 : theOutStatRec.define("median", outMdns);
1602 : }
1603 : }
1604 0 : return theOutStatRec;
1605 : }
1606 :
1607 : // robust image statistics for better noise estimation - revised for making it faster
1608 : // Modified version to do stats all at once
1609 6 : Record SDMaskHandler::calcRobustImageStatistics(ImageInterface<Float>& res, ImageInterface<Float>& prevmask , LatticeExpr<Bool>& pbmask, String& LELmask, Record* regionPtr, const Bool robust, Vector<Bool>& chanflag )
1610 : {
1611 18 : LogIO os( LogOrigin("SDMaskHandler","calcRobustImageStatistics",WHERE) );
1612 :
1613 : //Bool debugStats(true);
1614 :
1615 12 : Array<Float> wholemaskdata;
1616 12 : IPosition maskshp = prevmask.shape();
1617 12 : IPosition maskstart(4,0);
1618 12 : IPosition masklength(4,maskshp(0),maskshp(1), 1, 1);
1619 12 : Slicer masksl(maskstart, masklength);
1620 6 : prevmask.doGetSlice(wholemaskdata,masksl); // single plane
1621 6 : Float npixwholemask = sum(wholemaskdata);
1622 6 : Bool fullmask(False);
1623 :
1624 6 : if (npixwholemask == (Float) maskshp(0) * (Float) maskshp(1) ) {
1625 0 : fullmask=True;
1626 0 : os<<LogIO::DEBUG1 <<"Appears to be fully masked! npix for the whole mask ="<<npixwholemask<<LogIO::POST;
1627 : }
1628 :
1629 : //res.get(resdata);
1630 : //tempres->put(resdata);
1631 : // if input image (res) has a pixel mask, make sure to honor it so the region is exclude from statistics
1632 : //if (res.hasPixelMask()) {
1633 : // tempres->attachMask(res.pixelMask());
1634 : //}
1635 6 : if (nfalse(pbmask.getMask())) {
1636 0 : os<<LogIO::DEBUG1<<"has pbmask..."<<LogIO::POST;
1637 0 : os<<LogIO::DEBUG1<<"-> nfalse(pbmask)="<<nfalse(pbmask.getMask())<<LogIO::POST;
1638 : }
1639 :
1640 : // do stats on a whole cube at once for each algrithms
1641 : //check Stokes and spectral axes
1642 12 : IPosition shp = res.shape();
1643 6 : Int specaxis = CoordinateUtil::findSpectralAxis(res.coordinates());
1644 6 : uInt nchan = shp(specaxis);
1645 12 : Vector<Stokes::StokesTypes> whichPols;
1646 6 : Int stokesaxis = CoordinateUtil::findStokesAxis(whichPols, res.coordinates());
1647 6 : uInt nStokes= shp(stokesaxis);
1648 : uInt iaxis2;
1649 : uInt iaxis3;
1650 :
1651 : //TempImage<Bool> pbmaskim(shp, tempres->coordinates(), memoryToUse());
1652 : //TempImage<Bool> pbmaskim(shp, res.coordinates(), memoryToUse());
1653 : //pbmaskim.copyData(pbmask);
1654 :
1655 : //check dimensions of tempres, prevmask, pbmask
1656 6 : if (prevmask.shape()!=shp) {
1657 0 : throw(AipsError("Mismatch in shapes of the previous mask and residual image"));
1658 : }
1659 6 : else if (pbmask.shape()!=shp) {
1660 0 : throw(AipsError("Mismatch in shapes of pbmask and residual image"));
1661 : }
1662 :
1663 : //for stats storage
1664 12 : IPosition statshape(2, shp(2), shp(3));
1665 12 : Array<Double> outMins(statshape);
1666 12 : Array<Double> outMaxs(statshape);
1667 12 : Array<Double> outRmss(statshape);
1668 12 : Array<Double> outMads(statshape);
1669 12 : Array<Double> outMdns(statshape);
1670 :
1671 : //do stats for each algortims all at once.
1672 : // 2nd arg is regionRecord, 3rd is LELmask expression and those will be AND
1673 : // to define a region to be get statistics
1674 : //ImageStatsCalculator imcalc( tempres_ptr, 0, "", False);
1675 : // Chauvenet over a full image
1676 : // for an empty (= no mask) mask, use Chauvenet algorithm with maxiter=5 and zscore=-1
1677 6 : TempImage<Float>* tempRes = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
1678 6 : tempRes->copyData(res);
1679 6 : tempRes->attachMask(pbmask);
1680 12 : std::shared_ptr<casacore::ImageInterface<Float> > tempres_ptr(tempRes);
1681 12 : ImageStatsCalculator<Float> imcalc( tempres_ptr, regionPtr, LELmask, False);
1682 12 : Vector<Int> axes(2);
1683 6 : axes[0] = 0;
1684 6 : axes[1] = 1;
1685 6 : imcalc.setAxes(axes);
1686 6 : os<<LogIO::DEBUG1<<"Using Chauvenet algorithm for image statistics for a whole cube"<<LogIO::POST;
1687 6 : imcalc.configureChauvenet(Double(-1.0), Int(5));
1688 6 : imcalc.setRobust(robust);
1689 12 : Record thestatsNoMask = imcalc.statistics();
1690 :
1691 : // do stats outside the mask
1692 : //tempres_ptr.reset();
1693 12 : Record thestatsWithMask;
1694 6 : if(!fullmask) {
1695 6 : TempImage<Float>* tempRes2 = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
1696 6 : tempRes2->copyData(res);
1697 6 : os<<LogIO::DEBUG1<<"Do the image statitistics in a region outside the mask..."<<LogIO::POST;
1698 12 : LatticeExpr<Bool> outsideMaskReg;
1699 6 : outsideMaskReg = LatticeExpr<Bool> (iif(prevmask == 1.0 || !pbmask, False, True));
1700 : // need this case? for no existance of pb mask?
1701 : //outsideMaskReg = LatticeExpr<Bool> (iif(prevmask == 1.0, False, True));
1702 6 : tempRes2->attachMask(outsideMaskReg);
1703 12 : std::shared_ptr<casacore::ImageInterface<Float> > tempres_ptr2(tempRes2);
1704 6 : ImageStatsCalculator<Float> imcalc2( tempres_ptr2, regionPtr, LELmask, False);
1705 6 : imcalc2.setAxes(axes);
1706 6 : imcalc2.setRobust(robust);
1707 6 : thestatsWithMask = imcalc2.statistics();
1708 : }
1709 :
1710 12 : Array<Double> arrminsNoMask, arrmaxsNoMask, arrrmssNoMask, arrmadsNoMask, arrmdnsNoMask;
1711 12 : Array<Double> arrmins, arrmaxs, arrrmss, arrmads, arrmdns;
1712 6 : thestatsNoMask.get(RecordFieldId("min"), arrminsNoMask);
1713 6 : thestatsNoMask.get(RecordFieldId("max"), arrmaxsNoMask);
1714 6 : thestatsNoMask.get(RecordFieldId("rms"), arrrmssNoMask);
1715 12 : IPosition statsind(arrmins.ndim(), 0);
1716 : //robust = True only
1717 6 : if (robust) {
1718 6 : thestatsNoMask.get(RecordFieldId("medabsdevmed"), arrmadsNoMask);
1719 6 : thestatsNoMask.get(RecordFieldId("median"), arrmdnsNoMask);
1720 : }
1721 6 : if (!fullmask) { // with Mask stats
1722 6 : thestatsWithMask.get(RecordFieldId("min"), arrmins);
1723 6 : thestatsWithMask.get(RecordFieldId("max"), arrmaxs);
1724 6 : thestatsWithMask.get(RecordFieldId("rms"), arrrmss);
1725 6 : if (robust) {
1726 6 : thestatsWithMask.get(RecordFieldId("medabsdevmed"), arrmads);
1727 6 : thestatsWithMask.get(RecordFieldId("median"), arrmdns);
1728 : }
1729 : }
1730 :
1731 12 : for (uInt istokes = 0; istokes < nStokes; istokes++) {
1732 12 : std::vector<double> mins, maxs, rmss, mads, mdns;
1733 12 : for (uInt ichan = 0; ichan < nchan; ichan++ ) {
1734 6 : if (stokesaxis==3) {
1735 0 : iaxis2 = ichan;
1736 0 : iaxis3 = istokes;
1737 0 : statsind(0) = iaxis2;
1738 0 : if (nStokes>1) {
1739 0 : statsind(1) = iaxis3;
1740 : }
1741 : }
1742 : else {
1743 6 : iaxis2 = istokes;
1744 6 : iaxis3 = ichan;
1745 6 : if (nStokes>1) {
1746 0 : statsind(0) = iaxis2;
1747 0 : statsind(1) = iaxis3;
1748 : }
1749 : else { // Stokes axis is degenerate
1750 6 : statsind(0) = iaxis3;
1751 : }
1752 : }
1753 :
1754 6 : if (chanflag.nelements()==0 || !chanflag(ichan)) { // get new stats
1755 : // check if mask is empty (evaulated as the whole input mask )
1756 12 : Array<Float> maskdata;
1757 12 : IPosition start(4, 0, 0, iaxis2, iaxis3);
1758 12 : IPosition length(4, shp(0),shp(1), 1, 1);
1759 12 : Slicer sl(start, length);
1760 : //os<<"slicer for subimage start="<<start<<" length="<<length<<LogIO::POST;
1761 : // create subimage (slice) of tempres
1762 : //AxesSpecifier aspec(False);
1763 : //SubImage<Float>* subprevmask = new SubImage<Float>(prevmask, sl, aspec, True);
1764 :
1765 : // get chan mask data to evaulate no mask
1766 6 : prevmask.doGetSlice(maskdata,sl); // single plane
1767 6 : Float nmaskpix = sum(maskdata);
1768 12 : Array<Bool> booldata;
1769 :
1770 :
1771 : Double min, max, rms, mad, mdn;
1772 : // for an empty (= no mask) mask, use Chauvenet algorithm with maxiter=5 and zscore=-1
1773 6 : if (nmaskpix==0.0 || fullmask ) {
1774 4 : os<<"[C"<<ichan<<"] Using Chauvenet algorithm for the image statistics "<<LogIO::POST;
1775 4 : if(arrminsNoMask.nelements()==0 || arrmaxsNoMask.nelements() ==0 || arrrmssNoMask.nelements() ==0){
1776 0 : throw(AipsError("No image statistics possible on zero element data"));
1777 : }
1778 4 : min = arrminsNoMask(statsind);
1779 4 : max = arrmaxsNoMask(statsind);
1780 4 : rms = arrrmssNoMask(statsind);
1781 8 : if (robust) {
1782 4 : if(arrmadsNoMask.nelements()==0 || arrmdnsNoMask.nelements() ==0 || arrrmssNoMask.nelements() ==0)
1783 0 : throw(AipsError("No robust image statistics possible on zero element data"));
1784 4 : mad = arrmadsNoMask(statsind);
1785 4 : mdn = arrmdnsNoMask(statsind);
1786 : }
1787 : }
1788 : else {
1789 : //os<<"[C"<<ichan<<"] Using the image statitistics in a region outside the mask"<<LogIO::POST;
1790 2 : min = arrmins(statsind);
1791 2 : max = arrmaxs(statsind);
1792 2 : rms = arrrmss(statsind);
1793 2 : if (robust) {
1794 2 : mad = arrmads(statsind);
1795 2 : mdn = arrmdns(statsind);
1796 : }
1797 : }
1798 :
1799 :
1800 : // repack
1801 : //if (arrmaxs.nelements()==0 ){
1802 : // throw(AipsError("No image statisitics is returned. Possible the whole image is masked."));
1803 : //}
1804 :
1805 6 : mins.push_back(min);
1806 6 : maxs.push_back(max);
1807 6 : rmss.push_back(rms);
1808 6 : if (robust) {
1809 6 : mads.push_back(mad);
1810 6 : mdns.push_back(mdn);
1811 : }
1812 : }// if-ichanflag end
1813 : else {
1814 0 : mins.push_back(Double(0.0));
1815 0 : maxs.push_back(Double(0.0));
1816 0 : rmss.push_back(Double(0.0));
1817 0 : if (robust) {
1818 0 : mads.push_back(Double(0.0));
1819 0 : mdns.push_back(Double(0.0));
1820 : }
1821 : }
1822 : } // chan-for-loop end
1823 : //os<<" rms vector for stokes="<<istokes<<" : "<<rmss<<LogIO::POST;
1824 : //os<<"outMins.shape="<<outMins.shape()<<LogIO::POST;
1825 :
1826 12 : IPosition start(2, istokes, 0);
1827 12 : IPosition length(2, 1, nchan); // slicer per Stokes
1828 : //Slicer sl(blc, trc,Slicer::endIsLast );
1829 12 : Slicer slstokes(start, length);
1830 : //cerr<<"set outMin slstokes="<<slstokes<<" to the mins vector"<<endl;
1831 12 : Vector<Double> minvec(mins);
1832 12 : Vector<Double> maxvec(maxs);
1833 12 : Vector<Double> rmsvec(rmss);
1834 : //os<<"stl vector rmss"<<rmss<<LogIO::POST;
1835 : //os<<"rmsvec.shape="<<rmsvec.shape()<<" rmsvec="<<rmsvec<<LogIO::POST;
1836 12 : Matrix<Double> minmat(minvec);
1837 12 : Matrix<Double> maxmat(maxvec);
1838 : //Matrix<Double> rmsmat((Vector<Double>) rmss);
1839 12 : Matrix<Double> rmsmat(rmsvec);
1840 : //os<<"Intial shape of rmsmat="<<rmsmat.shape()<<LogIO::POST;
1841 : //os<<"Intial content of rmsmat="<<rmsmat<<LogIO::POST;
1842 : // copyValues do not work if there is a degerate axis in front
1843 6 : minmat.resize(IPosition(2, nchan, 1),True);
1844 6 : maxmat.resize(IPosition(2, nchan, 1),True);
1845 6 : rmsmat.resize(IPosition(2, nchan, 1),True);
1846 12 : Array<Double> minarr = transpose(minmat);
1847 12 : Array<Double> maxarr = transpose(maxmat);
1848 6 : Array<Double> rmsarr = transpose(rmsmat);
1849 6 : if (robust) {
1850 12 : Vector<Double> madvec(mads);
1851 12 : Vector<Double> mdnvec(mdns);
1852 12 : Matrix<Double> madmat(madvec);
1853 12 : Matrix<Double> mdnmat(mdnvec);
1854 6 : madmat.resize(IPosition(2, nchan, 1),True);
1855 6 : mdnmat.resize(IPosition(2, nchan, 1),True);
1856 12 : Array<Double> madarr = transpose(madmat);
1857 6 : Array<Double> mdnarr = transpose(mdnmat);
1858 6 : outMads(slstokes) = madarr;
1859 6 : outMdns(slstokes) = mdnarr;
1860 : //os<<"madarr="<<madarr<<LogIO::POST;
1861 : //os<<"mdnarr="<<mdnarr<<LogIO::POST;
1862 : }
1863 6 : outMins(slstokes) = minarr;
1864 6 : outMaxs(slstokes) = maxarr;
1865 6 : outRmss(slstokes) = rmsarr;
1866 : } // stokes-for-loop end
1867 :
1868 6 : Record theOutStatRec;
1869 6 : if (shp(2) == 1) { //Single Stokes plane
1870 : // remove degenerate axis
1871 : //os<<"Stokes axis has a single Stokes"<<LogIO::POST;
1872 6 : theOutStatRec.define("min", outMins.nonDegenerate(IPosition(1,1)));
1873 6 : theOutStatRec.define("max", outMaxs.nonDegenerate(IPosition(1,1)));
1874 6 : theOutStatRec.define("rms", outRmss.nonDegenerate(IPosition(1,1)));
1875 6 : if (robust) {
1876 6 : theOutStatRec.define("medabsdevmed", outMads.nonDegenerate(IPosition(1,1)));
1877 6 : theOutStatRec.define("median", outMdns.nonDegenerate(IPosition(1,1)));
1878 : }
1879 : }
1880 0 : else if(shp(3) == 1) { //Single channel plane
1881 0 : theOutStatRec.define("min", outMins.nonDegenerate(IPosition(1,0)));
1882 0 : theOutStatRec.define("max", outMaxs.nonDegenerate(IPosition(1,0)));
1883 0 : theOutStatRec.define("rms", outRmss.nonDegenerate(IPosition(1,0)));
1884 0 : if (robust) {
1885 0 : theOutStatRec.define("medabsdevmed", outMads.nonDegenerate(IPosition(1,0)));
1886 0 : theOutStatRec.define("median", outMdns.nonDegenerate(IPosition(1,0)));
1887 : }
1888 : }
1889 : else {
1890 0 : theOutStatRec.define("min", outMins);
1891 0 : theOutStatRec.define("max", outMaxs);
1892 0 : theOutStatRec.define("rms", outRmss);
1893 0 : if (robust) {
1894 0 : theOutStatRec.define("medabsdevmed", outMads);
1895 0 : theOutStatRec.define("median", outMdns);
1896 : }
1897 : }
1898 :
1899 12 : return theOutStatRec;
1900 : }
1901 :
1902 0 : void SDMaskHandler::autoMaskByThreshold(ImageInterface<Float>& mask,
1903 : const ImageInterface<Float>& res,
1904 : const ImageInterface<Float>& psf,
1905 : const Quantity& resolution,
1906 : const Float& resbybeam,
1907 : const Quantity& qthresh,
1908 : const Float& fracofpeak,
1909 : const Record& stats,
1910 : const Float& sigma,
1911 : const Int nmask,
1912 : const Bool autoadjust)
1913 : {
1914 0 : LogIO os( LogOrigin("SDMaskHandler","autoMaskByThreshold",WHERE) );
1915 0 : Array<Double> rms, max;
1916 : Double rmsthresh, minrmsval, maxrmsval, minmaxval, maxmaxval;
1917 0 : IPosition minrmspos, maxrmspos, minmaxpos, maxmaxpos;
1918 : Int npix;
1919 : //for debug set to True to save intermediate mask images on disk
1920 0 : Bool debug(False);
1921 :
1922 : //automask stage selecitons
1923 0 : Bool dobin(True);
1924 0 : Bool doregrid(True);
1925 0 : Bool doconvolve(True);
1926 :
1927 : // taking account for beam or input resolution
1928 0 : TempImage<Float> tempmask(mask.shape(), mask.coordinates(), memoryToUse());
1929 0 : IPosition shp = mask.shape();
1930 0 : CoordinateSystem incsys = res.coordinates();
1931 0 : Vector<Double> incVal = incsys.increment();
1932 0 : Vector<String> incUnit = incsys.worldAxisUnits();
1933 0 : Quantity qinc(incVal[0],incUnit[0]);
1934 0 : if (resolution.get().getValue() ) {
1935 : //npix = 2*Int(abs( resolution/(qinc.convert(resolution),qinc) ).getValue() );
1936 0 : npix = Int(abs( resolution/(qinc.convert(resolution),qinc) ).getValue() );
1937 0 : os << LogIO::NORMAL2 << "Use the input resolution:"<<resolution<<" fo binning "<< LogIO::POST;
1938 0 : os << LogIO::DEBUG1 << "inc = "<<qinc.getValue(resolution.getUnit())<<LogIO::POST;
1939 : }
1940 : else {
1941 : //use beam from residual or psf
1942 0 : ImageInfo resInfo = res.imageInfo();
1943 0 : ImageInfo psfInfo = psf.imageInfo();
1944 0 : GaussianBeam beam;
1945 0 : if (resInfo.hasBeam() || psfInfo.hasBeam()) {
1946 0 : if (resInfo.hasSingleBeam()) {
1947 0 : beam = resInfo.restoringBeam();
1948 : }
1949 0 : else if (resInfo.hasMultipleBeams()) {
1950 0 : beam = CasaImageBeamSet(resInfo.getBeamSet()).getCommonBeam();
1951 : }
1952 0 : else if (psfInfo.hasSingleBeam()) {
1953 0 : beam = psfInfo.restoringBeam();
1954 : }
1955 : else {
1956 0 : beam = CasaImageBeamSet(psfInfo.getBeamSet()).getCommonBeam();
1957 : }
1958 0 : Quantity bmaj = beam.getMajor();
1959 0 : Quantity bmin = beam.getMinor();
1960 0 : if (resbybeam > 0.0 ) {
1961 : //npix = 2*Int( Double(resbybeam) * abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
1962 0 : npix = Int( Double(resbybeam) * abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
1963 0 : os << LogIO::NORMAL2 << "Use "<< resbybeam <<" x beam size(maj)="<< Double(resbybeam)*bmaj <<" for binning."<< LogIO::POST;
1964 : }
1965 : else {
1966 : //npix = 2*Int( abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
1967 0 : npix = Int( abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
1968 0 : os << LogIO::NORMAL2 << "Use a beam size(maj):"<<bmaj<<" for binning."<< LogIO::POST;
1969 : }
1970 : }
1971 : else {
1972 0 : throw(AipsError("No restoring beam(s) in the input image/psf or resolution is given"));
1973 : }
1974 : }
1975 0 : os << LogIO::DEBUG1 << "Acutal bin size used: npix="<<npix<< LogIO::POST;
1976 0 : if (npix==0) {
1977 0 : os << "Resolution too small. No binning (nbin=1) is applied input image to evaluate the threshold." << LogIO::POST;
1978 0 : npix=1;
1979 : }
1980 :
1981 : // Determine threshold from input image stats
1982 0 : stats.get(RecordFieldId("max"), max);
1983 0 : stats.get(RecordFieldId("rms"), rms);
1984 0 : minMax(minmaxval,maxmaxval,minmaxpos, maxmaxpos, max);
1985 0 : minMax(minrmsval,maxrmsval,minrmspos, maxrmspos, rms);
1986 0 : os << LogIO::DEBUG1 <<"stats on the image: max="<<maxmaxval<<" rms="<<maxrmsval<<endl;
1987 0 : if (fracofpeak) {
1988 0 : rmsthresh = maxmaxval * fracofpeak;
1989 : //os << LogIO::NORMAL <<"Threshold by fraction of the peak(="<<fracofpeak<<") * max: "<<rmsthresh<< LogIO::POST;
1990 0 : os << LogIO::DEBUG1 <<"max at "<<maxmaxpos<<", dynamic range = "<<maxmaxval/rms(maxmaxpos) << LogIO::POST;
1991 : }
1992 0 : else if (sigma) {
1993 : //cerr<<"minval="<<minval<<" maxval="<<maxval<<endl;
1994 0 : rmsthresh = maxrmsval * sigma;
1995 : //os << LogIO::NORMAL <<"Threshold by sigma(="<<sigma<<")* rms (="<<maxrmsval<<") :"<<rmsthresh<< LogIO::POST;
1996 0 : os << LogIO::DEBUG1 <<"max rms at "<<maxrmspos<<", dynamic range = "<<max(maxrmspos)/maxrmsval << LogIO::POST;
1997 : }
1998 : else {
1999 0 : rmsthresh = qthresh.getValue(Unit("Jy"));
2000 0 : if ( rmsthresh==0.0 )
2001 0 : { throw(AipsError("Threshold for automask is not set"));}
2002 : }
2003 : //os << LogIO::NORMAL2 <<" thresh="<<rmsthresh<<LogIO::POST;
2004 :
2005 :
2006 0 : TempImage<Float>* tempIm2 = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
2007 0 : TempImage<Float>* tempIm = new TempImage<Float>(res.shape(), res.coordinates(), memoryToUse());
2008 0 : tempIm->copyData(res);
2009 :
2010 0 : SPCIIF tempIm2_ptr(tempIm2);
2011 0 : SPIIF tempIm3_ptr(tempIm);
2012 0 : SPIIF tempIm_ptr;
2013 : //
2014 : //binning stage
2015 0 : if (dobin) {
2016 0 : tempIm_ptr = makeMaskFromBinnedImage(res, npix, npix, fracofpeak, sigma, nmask, autoadjust, rmsthresh);
2017 : //for debugging: save the mask at this stage
2018 0 : if (debug) {
2019 0 : PagedImage<Float> tempBinIm(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(), "binnedThresh.Im");
2020 0 : tempBinIm.copyData(*(tempIm_ptr.get()));
2021 : }
2022 : }
2023 0 : if (doregrid) {
2024 : //regrid
2025 0 : os << LogIO::DEBUG1 <<" now regridding..."<<LogIO::POST;
2026 0 : IPosition axes(3,0, 1, 2);
2027 0 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(incsys);
2028 0 : axes(0) = dirAxes(0);
2029 0 : axes(1) = dirAxes(1);
2030 0 : axes(2) = CoordinateUtil::findSpectralAxis(incsys);
2031 0 : Record* dummyrec = 0;
2032 0 : SPCIIF inmask_ptr(tempIm_ptr);
2033 0 : ImageRegridder<Float> regridder(inmask_ptr, "", tempIm2_ptr, axes, dummyrec, "", True, shp);
2034 0 : regridder.setMethod(Interpolate2D::LINEAR);
2035 0 : tempIm_ptr = regridder.regrid();
2036 : //for debugging: save the mask at this stage
2037 0 : if (debug) {
2038 0 : PagedImage<Float> tempGridded(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(), "binAndGridded.Im");
2039 0 : tempGridded.copyData(*(tempIm_ptr.get()));
2040 : }
2041 : }
2042 : else {
2043 0 : tempIm_ptr = tempIm3_ptr;
2044 : }
2045 0 : if (doconvolve) {
2046 : //
2047 0 : SPIIF outmask = convolveMask(*(tempIm_ptr.get()), npix, npix);
2048 0 : tempIm_ptr = outmask;
2049 : //
2050 : //for debugging: save the mask at this stage
2051 0 : if (debug) {
2052 0 : PagedImage<Float> tempconvIm(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(),"convolved.Im");
2053 0 : tempconvIm.copyData(*(tempIm_ptr.get()));
2054 : }
2055 : //os<<"done convolving the mask "<<LogIO::POST;
2056 : }
2057 :
2058 : //os <<"Final thresholding with rmsthresh/afactor="<< rmsthresh/afactor <<LogIO::POST;
2059 : //LatticeExpr<Float> themask( iif( *(tempIm_ptr.get()) > rmsthresh/afactor, 1.0, 0.0 ));
2060 : // previous 1/0 mask (regridded), max pix value should be <1.0, take arbitary cut off at 0.1
2061 0 : LatticeExpr<Float> themask( iif( *(tempIm_ptr.get()) > 0.1, 1.0, 0.0 ));
2062 0 : if (res.hasPixelMask()) {
2063 0 : LatticeExpr<Bool> pixmask(res.pixelMask());
2064 0 : mask.copyData( (LatticeExpr<Float>)( iif((mask + themask) > 0.0 && pixmask, 1.0, 0.0 ) ) );
2065 0 : os <<LogIO::DEBUG1 <<"add previous mask, pbmask and the new mask.."<<LogIO::POST;
2066 : }
2067 : else {
2068 : //for debug
2069 : /***
2070 : PagedImage<Float> tempthemask(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(),"tempthemask.Im");
2071 : tempthemask.copyData(themask);
2072 : ***/
2073 :
2074 : //os <<"Lattice themask is created..."<<LogIO::POST;
2075 : //LatticeExpr<Float> themask( iif( tempconvim > rmsthresh/afactor, 1.0, 0.0 ));
2076 0 : mask.copyData( (LatticeExpr<Float>)( iif((mask + themask) > 0.0, 1.0, 0.0 ) ) );
2077 0 : os <<LogIO::DEBUG1 <<"add previous mask and the new mask.."<<LogIO::POST;
2078 : }
2079 0 : }
2080 :
2081 0 : void SDMaskHandler::autoMaskByThreshold2(ImageInterface<Float>& mask,
2082 : const ImageInterface<Float>& res,
2083 : const ImageInterface<Float>& psf,
2084 : const Quantity& resolution,
2085 : const Float& resbybeam,
2086 : const Quantity& qthresh,
2087 : const Float& fracofpeak,
2088 : const Record& stats,
2089 : const Float& sigma,
2090 : const Int nmask)
2091 : {
2092 0 : LogIO os( LogOrigin("SDMaskHandler","autoMaskByThreshold2",WHERE) );
2093 0 : Array<Double> rms, max;
2094 : Double rmsthresh, minrmsval, maxrmsval, minmaxval, maxmaxval;
2095 0 : IPosition minrmspos, maxrmspos, minmaxpos, maxmaxpos;
2096 : Int npix;
2097 : Int beampix;
2098 :
2099 : //for debug set to True to save intermediate mask images on disk
2100 : //Bool debug(False);
2101 :
2102 : // taking account for beam or input resolution
2103 0 : TempImage<Float> tempmask(mask.shape(), mask.coordinates(), memoryToUse());
2104 0 : IPosition shp = mask.shape();
2105 0 : CoordinateSystem incsys = res.coordinates();
2106 0 : Vector<Double> incVal = incsys.increment();
2107 0 : Vector<String> incUnit = incsys.worldAxisUnits();
2108 0 : Quantity qinc(incVal[0],incUnit[0]);
2109 0 : if (resolution.get().getValue() ) {
2110 0 : npix = Int(abs( resolution/(qinc.convert(resolution),qinc) ).getValue() );
2111 0 : beampix = Int( C::pi * npix * npix /(4.*C::ln2));
2112 0 : os << LogIO::NORMAL2 << "Use the input resolution:"<<resolution<<" for pruning "<< LogIO::POST;
2113 0 : os << LogIO::DEBUG1 << "inc = "<<qinc.getValue(resolution.getUnit())<<LogIO::POST;
2114 : }
2115 : else {
2116 : //use beam from residual or psf
2117 0 : ImageInfo resInfo = res.imageInfo();
2118 0 : ImageInfo psfInfo = psf.imageInfo();
2119 0 : GaussianBeam beam;
2120 : Int npixmin;
2121 0 : if (resInfo.hasBeam() || psfInfo.hasBeam()) {
2122 0 : if (resInfo.hasSingleBeam()) {
2123 0 : beam = resInfo.restoringBeam();
2124 : }
2125 0 : else if (resInfo.hasMultipleBeams()) {
2126 0 : beam = CasaImageBeamSet(resInfo.getBeamSet()).getCommonBeam();
2127 : }
2128 0 : else if (psfInfo.hasSingleBeam()) {
2129 0 : beam = psfInfo.restoringBeam();
2130 : }
2131 : else {
2132 0 : beam = CasaImageBeamSet(psfInfo.getBeamSet()).getCommonBeam();
2133 : }
2134 0 : Quantity bmaj = beam.getMajor();
2135 0 : Quantity bmin = beam.getMinor();
2136 0 : if (resbybeam > 0.0 ) {
2137 0 : npix = Int( Double(resbybeam) * abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
2138 0 : npixmin = Int( Double(resbybeam) * abs( (bmin/(qinc.convert(bmin),qinc)).get().getValue() ) );
2139 0 : beampix = Int(C::pi * npix * npixmin / (4. * C::ln2));
2140 :
2141 0 : os << LogIO::NORMAL2 << "Use "<< resbybeam <<" x beam size(maj)="<< Double(resbybeam)*bmaj <<" for pruning."<< LogIO::POST;
2142 : }
2143 : else {
2144 0 : npix = Int( abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
2145 0 : npixmin = Int( abs( (bmin/(qinc.convert(bmin),qinc)).get().getValue() ) );
2146 0 : beampix = Int(C::pi * npix * npixmin / (4. * C::ln2));
2147 0 : os << LogIO::NORMAL2 << "Use a beam size(maj):"<<bmaj<<" for pruning."<< LogIO::POST;
2148 : }
2149 : }
2150 : else {
2151 0 : throw(AipsError("No restoring beam(s) in the input image/psf or resolution is given"));
2152 : }
2153 : }
2154 0 : os << LogIO::DEBUG1 << "Acutal bin size used: npix="<<npix<< LogIO::POST;
2155 :
2156 : // Determine threshold from input image stats
2157 0 : stats.get(RecordFieldId("max"), max);
2158 0 : stats.get(RecordFieldId("rms"), rms);
2159 0 : minMax(minmaxval,maxmaxval,minmaxpos, maxmaxpos, max);
2160 0 : minMax(minrmsval,maxrmsval,minrmspos, maxrmspos, rms);
2161 0 : os << LogIO::DEBUG1 <<"stats on the image: max="<<maxmaxval<<" rms="<<maxrmsval<<endl;
2162 0 : if (fracofpeak) {
2163 0 : rmsthresh = maxmaxval * fracofpeak;
2164 0 : os << LogIO::NORMAL <<"Threshold by fraction of the peak(="<<fracofpeak<<") * max: "<<rmsthresh<< LogIO::POST;
2165 0 : os << LogIO::DEBUG1 <<"max at "<<maxmaxpos<<", dynamic range = "<<maxmaxval/rms(maxmaxpos) << LogIO::POST;
2166 : }
2167 0 : else if (sigma) {
2168 : //cerr<<"minval="<<minval<<" maxval="<<maxval<<endl;
2169 0 : rmsthresh = maxrmsval * sigma;
2170 0 : os << LogIO::NORMAL <<"Threshold by sigma(="<<sigma<<")* rms (="<<maxrmsval<<") :"<<rmsthresh<< LogIO::POST;
2171 0 : os << LogIO::DEBUG1 <<"max rms at "<<maxrmspos<<", dynamic range = "<<max(maxrmspos)/maxrmsval << LogIO::POST;
2172 : }
2173 : else {
2174 0 : rmsthresh = qthresh.getValue(Unit("Jy"));
2175 0 : if ( rmsthresh==0.0 )
2176 0 : { throw(AipsError("Threshold for automask is not set"));}
2177 : }
2178 0 : os << LogIO::NORMAL2 <<" thresh="<<rmsthresh<<LogIO::POST;
2179 :
2180 0 : std::shared_ptr<ImageInterface<Float> > tempIm_ptr = pruneRegions(res, rmsthresh, nmask, beampix);
2181 0 : LatticeExpr<Float> themask( iif( *(tempIm_ptr.get()) > rmsthresh, 1.0, 0.0 ));
2182 :
2183 : //for debug
2184 : /***
2185 : PagedImage<Float> tempthemask(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(),"tempthemask.Im");
2186 : tempthemask.copyData(themask);
2187 : ***/
2188 0 : if (res.hasPixelMask()) {
2189 0 : LatticeExpr<Bool> pixmask(res.pixelMask());
2190 0 : mask.copyData( (LatticeExpr<Float>)( iif((mask + themask) > 0.0 && pixmask, 1.0, 0.0 ) ) );
2191 0 : mask.clearCache();
2192 0 : mask.unlock();
2193 0 : mask.tempClose();
2194 0 : os <<LogIO::DEBUG1 <<"Add previous mask, pbmask and the new mask.."<<LogIO::POST;
2195 : }
2196 : else {
2197 : //os <<"Lattice themask is created..."<<LogIO::POST;
2198 : //LatticeExpr<Float> themask( iif( tempconvim > rmsthresh/afactor, 1.0, 0.0 ));
2199 0 : mask.copyData( (LatticeExpr<Float>)( iif((mask + themask) > 0.0, 1.0, 0.0 ) ) );
2200 0 : os <<LogIO::DEBUG1 <<"Add previous mask and the new mask.."<<LogIO::POST;
2201 : }
2202 0 : }//end of makeAutoMaskByThreshold2
2203 :
2204 : // *** auto-multithresh ***
2205 : // for implemtation of Amanda's algorithm
2206 52 : void SDMaskHandler::autoMaskByMultiThreshold(ImageInterface<Float>& mask,
2207 : ImageInterface<Float>& posmask,
2208 : const ImageInterface<Float>& res,
2209 : const ImageInterface<Float>& psf,
2210 : const Record& stats,
2211 : const Record& robuststats,
2212 : const Int iterdone,
2213 : Vector<Bool>& chanFlag,
2214 : const Float& minPercentChange,
2215 : const Float& sidelobeLevel,
2216 : const Float& sidelobeThresholdFactor,
2217 : const Float& noiseThresholdFactor,
2218 : const Float& lowNoiseThresholdFactor,
2219 : const Float& negativeThresholdFactor,
2220 : const Float& cutThreshold,
2221 : const Float& smoothFactor,
2222 : const Float& minBeamFrac,
2223 : const Int growIterations,
2224 : const Bool doGrowPrune,
2225 : const Bool verbose,
2226 : const Bool isthresholdreached)
2227 : {
2228 156 : LogIO os( LogOrigin("SDMaskHandler","autoMaskByMultiThreshold",WHERE) );
2229 104 : Array<Double> rmss, maxs, mins, mads, mdns;
2230 104 : Array<Double> resRmss;
2231 104 : IPosition minrmspos, maxrmspos, minmaxpos, maxmaxpos, minmadpos, maxmadpos;
2232 : Int nxpix, nypix;
2233 :
2234 : // % min mask pixel change (to trigger new automask creation per chan plane) to a fractional change
2235 52 : Float fracChange = minPercentChange/100.0;
2236 :
2237 : //store summary info
2238 104 : Record summaryRec;
2239 52 : summaryRec.define("sidelobelevel",sidelobeLevel);
2240 :
2241 : //for debug set to True to save intermediate mask images on disk
2242 52 : Bool debug(false); // create additional temp masks for debugging
2243 52 : Bool debug2(false); // debug2 saves masks before/after prune and binary dilation
2244 :
2245 : //set true to use calcImageStatistics2 and thresholds adjusted for the location (median)
2246 : //Bool newstats(True); // turn on new stats definition of threshold calc.
2247 :
2248 : //Timer
2249 52 : Timer timer;
2250 :
2251 : //debug
2252 52 : if (debug2) {
2253 0 : PagedImage<Float> tempcurinmask(mask.shape(), mask.coordinates(), "currrent-in-mask-"+String::toString(iterdone)+".im");
2254 0 : tempcurinmask.copyData(mask);
2255 : }
2256 :
2257 : // tempmsk: working image for the curret mask
2258 156 : TempImage<Float> tempmask(mask.shape(), mask.coordinates(), memoryToUse());
2259 52 : tempmask.set(0);
2260 : // prevmask: mask from previous iter.
2261 156 : TempImage<Float> prevmask(mask.shape(), mask.coordinates(), memoryToUse());
2262 : // use positive only previous mask
2263 : //prevmask.copyData(LatticeExpr<Float>(mask) );
2264 52 : prevmask.copyData(LatticeExpr<Float>(posmask) );
2265 : // set up a container for a full cube negative mask
2266 : //if (negativeThresholdFactor > 0) {
2267 156 : TempImage<Float> thenegmask(mask.shape(), mask.coordinates(), memoryToUse());
2268 52 : thenegmask.set(0);
2269 : //}
2270 : // taking account for beam or input resolution
2271 : //IPosition shp = mask.shape();
2272 104 : CoordinateSystem incsys = res.coordinates();
2273 104 : Vector<Double> incVal = incsys.increment();
2274 104 : Vector<String> incUnit = incsys.worldAxisUnits();
2275 104 : Quantity qinc(incVal[0],incUnit[0]);
2276 : //use beam from residual or psf
2277 104 : ImageInfo resInfo = res.imageInfo();
2278 104 : ImageInfo psfInfo = psf.imageInfo();
2279 :
2280 104 : GaussianBeam beam, modbeam; // modbeam for smooth
2281 : Double pruneSize;
2282 52 : if (resInfo.hasBeam() || psfInfo.hasBeam()) {
2283 52 : if (resInfo.hasSingleBeam()) {
2284 0 : beam = resInfo.restoringBeam();
2285 : }
2286 52 : else if (resInfo.hasMultipleBeams()) {
2287 0 : beam = CasaImageBeamSet(resInfo.getBeamSet()).getCommonBeam();
2288 : }
2289 52 : else if (psfInfo.hasSingleBeam()) {
2290 40 : beam = psfInfo.restoringBeam();
2291 : }
2292 : else {
2293 12 : beam = CasaImageBeamSet(psfInfo.getBeamSet()).getCommonBeam();
2294 : }
2295 : // check
2296 52 : if(std::isinf( beam.getMajor(Unit("arcsec"))) || std::isinf( beam.getMinor(Unit("arcsec"))) ){
2297 0 : throw(AipsError("A bad common beam, which is used to set smoothing and pruning sizes for automask. At least one of the axes of the beam is infinite."));
2298 : }
2299 104 : Quantity bmaj = beam.getMajor();
2300 52 : Quantity bmin = beam.getMinor();
2301 : //for pruning for now
2302 : // minBeamFrac * beamarea
2303 52 : Double beamfrac=1.0;
2304 52 : if (minBeamFrac > 0.0) {
2305 52 : beamfrac = (Double) minBeamFrac;
2306 : }
2307 52 : Double beampix = pixelBeamArea(beam, incsys);
2308 52 : pruneSize = beamfrac * beampix;
2309 : //beam in pixels
2310 52 : nxpix = Int( Double(smoothFactor) * abs( (bmaj/(qinc.convert(bmaj),qinc)).get().getValue() ) );
2311 52 : nypix = Int( Double(smoothFactor) * abs( (bmin/(qinc.convert(bmin),qinc)).get().getValue() ) );
2312 52 : modbeam.setMajorMinor(Double(smoothFactor) * bmaj, Double(smoothFactor) * bmin);
2313 52 : modbeam.setPA(beam.getPA());
2314 :
2315 52 : os<<LogIO::DEBUG1<<"beam in pixels: B_maj="<<nxpix<<" B_min="<<nypix<<" beam area="<<beampix<<LogIO::POST;
2316 52 : os<<LogIO::NORMAL <<"prune size="<<pruneSize<<"(minbeamfrac="<<minBeamFrac<<" * beampix="<<beampix<<")"<<LogIO::POST;
2317 52 : summaryRec.define("pruneregionsize",pruneSize);
2318 : }
2319 : else {
2320 0 : throw(AipsError("No restoring beam(s) in the input image/psf"));
2321 : }
2322 :
2323 :
2324 : //One time parameter checks
2325 52 : if (!iterdone) {
2326 31 : if (cutThreshold >=0.2) {
2327 1 : os<<LogIO::WARN<<"Faint regions may not be included in the final mask. Consider decreasing cutthreshold."<<LogIO::POST;
2328 : }
2329 : }
2330 :
2331 : // Determine threshold from input image stats
2332 52 : stats.get(RecordFieldId("max"), maxs);
2333 52 : stats.get(RecordFieldId("min"), mins);
2334 : // robuststats may contains only robustrms and medians if it is filled by calcRobustRMS
2335 52 : if (robuststats.isDefined("rms")) {
2336 48 : robuststats.get(RecordFieldId("rms"), rmss);
2337 48 : robuststats.get(RecordFieldId("medabsdevmed"), mads);
2338 48 : resRmss = mads * 1.4826;
2339 48 : os<<LogIO::DEBUG1<<" rms from MAD (mads*1.4826)= "<<resRmss<<LogIO::POST;
2340 : }
2341 4 : else if (robuststats.isDefined("robustrms") ) {
2342 4 : robuststats.get(RecordFieldId("robustrms"),resRmss); // already converted from MAD to rms
2343 4 : os<<LogIO::DEBUG1<<" robustrms from MAD (mads*1.4826)= "<<resRmss<<LogIO::POST;
2344 : }
2345 52 : os<<LogIO::DEBUG1<<"get mdns"<<LogIO::POST;
2346 52 : robuststats.get(RecordFieldId("median"), mdns);
2347 :
2348 : //check for pbmask
2349 104 : IPosition imshp=res.shape();
2350 104 : IPosition imstart(4, 0, 0, 0, 0);
2351 104 : IPosition imlength(4, imshp(0),imshp(1), imshp(2), imshp(3));
2352 :
2353 : Float sidelobeThreshold;
2354 : Float noiseThreshold;
2355 : Float lowNoiseThreshold;
2356 : Float negativeThreshold;
2357 : //
2358 : //determine shape, nchan, npol from residual image
2359 104 : CoordinateSystem imcoord = res.coordinates();
2360 52 : Int specAxis = CoordinateUtil::findSpectralAxis(imcoord);
2361 104 : Vector <Stokes::StokesTypes> whichPols;
2362 52 : Int polAxis = CoordinateUtil::findStokesAxis(whichPols, imcoord);
2363 52 : Int nchan = -1;
2364 52 : Int npol = -1;
2365 52 : if (specAxis != -1) {
2366 52 : nchan = imshp(specAxis);
2367 : }
2368 52 : if (polAxis != -1 ) {
2369 52 : npol = imshp(polAxis);
2370 : }
2371 : //Int specAxis = CoordinateUtil::findSpectralAxis(res.coordinates());
2372 : // here, now chindx really means index to extract per-plane
2373 : //
2374 104 : IPosition statshp = mdns.shape();
2375 104 : IPosition chindx = statshp;
2376 52 : Int poldim = (npol == -1? 1:npol);
2377 52 : Int chandim = (nchan == -1? 1:nchan);
2378 104 : Matrix<Float> maskThreshold(poldim, chandim);
2379 104 : Matrix<Float> lowMaskThreshold(poldim, chandim);
2380 104 : Matrix<Float> negativeMaskThreshold(poldim, chandim);
2381 104 : Matrix<String> ThresholdType(poldim, chandim);
2382 104 : Matrix<Bool> pruned(poldim, chandim);
2383 :
2384 : //for (uInt ich=0; ich < mads.nelements(); ich++) {
2385 260 : for (uInt ich=0; ich < (uInt)nchan; ich++) {
2386 : // for loop for Stokes as well?
2387 440 : for (uInt ipol=0; ipol < (uInt)npol; ipol++) {
2388 232 : if (nchan!=-1) {
2389 232 : if(npol==-1 || npol==1) {
2390 200 : chindx(0) = ich;
2391 : }
2392 : else {
2393 32 : chindx(0) = ipol;
2394 32 : chindx(1) = ich;
2395 : }
2396 : }
2397 : else { // pol only
2398 0 : chindx(0) = ipol;
2399 : }
2400 :
2401 : // turn on a new definition for new stats --- remove old one once tested
2402 : //if (newstats) {
2403 : // os<<LogIO::DEBUG1<<"Using the new statistics ..."<<LogIO::POST;
2404 : // sidelobeThreshold = (Float)mdns(chindx) + sidelobeLevel * sidelobeThresholdFactor * (Float)maxs(chindx);
2405 : //}
2406 : //else {
2407 : // sidelobeThreshold = sidelobeLevel * sidelobeThresholdFactor * (Float)maxs(chindx);
2408 : //}
2409 :
2410 : // turn on a new definition for new stats --- remove old one once tested
2411 : //if (newstats) {
2412 : // noiseThreshold = (Float)mdns(chindx) + noiseThresholdFactor * (Float)resRmss(chindx);
2413 : // lowNoiseThreshold = (Float)mdns(chindx) + lowNoiseThresholdFactor * (Float)resRmss(chindx);
2414 : // negativeThreshold = (Float)mdns(chindx) + negativeThresholdFactor * (Float)resRmss(chindx);
2415 : //}
2416 : //else {
2417 : // noiseThreshold = noiseThresholdFactor * (Float)resRmss(chindx);
2418 : // lowNoiseThreshold = lowNoiseThresholdFactor * (Float)resRmss(chindx);
2419 : // negativeThreshold = negativeThresholdFactor * (Float)resRmss(chindx);
2420 : //}
2421 :
2422 : // include location (=median) for both fastnoise=true and false
2423 232 : Double absmax = max(abs(maxs(chindx)), abs(mins(chindx)));
2424 : // start with no offset
2425 232 : sidelobeThreshold = sidelobeLevel * sidelobeThresholdFactor * (Float)absmax;
2426 232 : noiseThreshold = noiseThresholdFactor * (Float)resRmss(chindx);
2427 232 : lowNoiseThreshold = lowNoiseThresholdFactor * (Float)resRmss(chindx);
2428 232 : negativeThreshold = negativeThresholdFactor * (Float)resRmss(chindx);
2429 : //negativeMaskThreshold(ich) = (-1.0)*max(sidelobeThreshold, negativeThreshold) + (Float)mdns(chindx);
2430 232 : negativeMaskThreshold(ipol, ich) = (-1.0)*max(sidelobeThreshold, negativeThreshold) + (Float)mdns(chindx);
2431 : // add the offset
2432 232 : sidelobeThreshold += (Float)mdns(chindx);
2433 232 : noiseThreshold += (Float)mdns(chindx);
2434 232 : lowNoiseThreshold += (Float)mdns(chindx);
2435 232 : negativeThreshold += (Float)mdns(chindx);
2436 : //maskThreshold(ich) = max(sidelobeThreshold, noiseThreshold);
2437 232 : maskThreshold(ipol, ich) = max(sidelobeThreshold, noiseThreshold);
2438 : //lowMaskThreshold(ich) = max(sidelobeThreshold, lowNoiseThreshold);
2439 232 : lowMaskThreshold(ipol, ich) = max(sidelobeThreshold, lowNoiseThreshold);
2440 : //ThresholdType(ich) = (maskThreshold(ich) == sidelobeThreshold? "sidelobe": "noise");
2441 232 : ThresholdType(ipol, ich) = (maskThreshold(ipol, ich) == sidelobeThreshold? "sidelobe": "noise");
2442 :
2443 232 : os << LogIO::DEBUG1 <<" sidelobeTreshold="<<sidelobeThreshold<<" noiseThreshold="<<noiseThreshold<<" lowNoiseThreshold="<<lowNoiseThreshold<<LogIO::POST;
2444 232 : os << LogIO::DEBUG1 <<" negativeThreshold(abs)="<<negativeThreshold<<", all thresholds include location ="<<(Float)mdns(chindx)<<LogIO::POST;
2445 : //os << LogIO::DEBUG1 <<" Using "<<ThresholdType(ich)<<" threshold for chan "<<String::toString(ich)<<" threshold="<<maskThreshold(ich)<<LogIO::POST;
2446 232 : os << LogIO::DEBUG1 <<" Using "<<ThresholdType(ipol, ich)<<" threshold for pol "<<String::toString(ipol)<<", chan "<< String::toString(ich)<<" threshold="<<maskThreshold(ipol, ich)<<LogIO::POST;
2447 : } // for-ipol
2448 : } // for-ich
2449 :
2450 :
2451 : //main per plane loop start here
2452 104 : IPosition planeshp(imshp.nelements(), 1);
2453 52 : planeshp(0) = imshp(0);
2454 52 : planeshp(1) = imshp(1);
2455 :
2456 : // store in matrix instead of vector to support full pols
2457 104 : Matrix<uInt> nreg(poldim, chandim,0);
2458 104 : Matrix<uInt> npruned(poldim, chandim, 0);
2459 104 : Vector<Float> dummysizes;
2460 104 : Matrix<uInt> ngrowreg(poldim, chandim, 0);
2461 104 : Matrix<uInt> ngrowpruned(poldim, chandim, 0);
2462 104 : Matrix<Float> negmaskpixs(poldim, chandim ,0);
2463 104 : Matrix<Bool> allPruned(poldim, chandim, False);
2464 : // for timing info storage
2465 104 : Vector<Double> timeInitThresh(3,0.0);
2466 104 : Vector<Double> timePrune(3,0.0);
2467 104 : Vector<Double> timeSmooth(3,0.0);
2468 104 : Vector<Double> timeGrow(3,0.0);
2469 104 : Vector<Double> timePruneGrow(3,0.0);
2470 104 : Vector<Double> timeSmoothGrow(3,0.0);
2471 104 : Vector<Double> timeNegThresh(3,0.0);
2472 :
2473 52 : Bool perplanetiming(True);
2474 :
2475 260 : for (uInt ich=0; ich < (uInt)nchan; ich++) {
2476 208 : if (npol==1) {
2477 200 : os << LogIO::NORMAL<< "*** Start auto-multithresh processing for Channel "<<ich<<"***"<<LogIO::POST;
2478 : }
2479 440 : for (uInt ipol=0; ipol < (uInt)npol; ipol++ ) {
2480 232 : if (npol!=1) {
2481 32 : os << LogIO::NORMAL<< "*** Start auto-multithresh processing for Channel "<<ich<<", Polarization "<<ipol<<"***"<<LogIO::POST;
2482 : }
2483 : // channel skip check
2484 232 : if (chanFlag(ich)) {
2485 0 : os << LogIO::NORMAL<<" Skip this channel "<<LogIO::POST;
2486 : }
2487 : else {
2488 : // Below corresponds to createThresholdMask in Amanda's Python code.
2489 464 : IPosition start(planeshp.nelements(),0);
2490 232 : if (specAxis != -1) {
2491 232 : start(specAxis)=ich;
2492 : }
2493 232 : if (polAxis != -1) {
2494 232 : start(polAxis)=ipol;
2495 : }
2496 :
2497 464 : IPosition length(planeshp.nelements(), planeshp(0), planeshp(1), 1, 1);
2498 464 : Slicer sl(start, length);
2499 464 : AxesSpecifier aspec(True); // keep degenerate axes
2500 464 : SubImage<Float> planeResImage(res, sl, aspec, true);
2501 696 : TempImage<Float> planeTempMask(planeResImage.shape(), planeResImage.coordinates(), memoryToUse());
2502 232 : planeTempMask.set(0); // initialize
2503 : //SubImage<Float> subprevmask(prevmask, sl, true, aspec, true);
2504 : // working copy of per-plane previous total mask to be modified. Started with en empty mask.
2505 : // For an un-touched version of per-plane previous mask, subMask is created at the grow mask step.
2506 464 : TempImage<Float> subprevmask(planeshp, planeResImage.coordinates(), memoryToUse());
2507 232 : subprevmask.set(0);
2508 464 : SubImage<Float> subposmask(posmask, sl, true, aspec, true);
2509 : //Vector<Bool> allPruned(nchan);
2510 : // sigle element vectors for input
2511 464 : Vector<Bool> chanFlag1elem(1);
2512 232 : chanFlag1elem(0) = chanFlag(ich);
2513 464 : Vector<Float> maskThreshold1elem(1);
2514 232 : maskThreshold1elem(0) = maskThreshold(ipol,ich);
2515 464 : Vector<Float> lowMaskThreshold1elem(1);
2516 232 : lowMaskThreshold1elem(0) = lowMaskThreshold(ipol,ich);
2517 464 : Vector<Float> negativeMaskThreshold1elem(1);
2518 232 : negativeMaskThreshold1elem(0) = negativeMaskThreshold(ipol,ich);
2519 :
2520 : // *** Pruning ***
2521 232 : if (minBeamFrac > 0.0 ) {
2522 : // do pruning...
2523 : //os<<LogIO::NORMAL<<"Pruning the current mask"<<LogIO::POST;
2524 232 : os << LogIO::NORMAL3 << "Start thresholding: create an initial mask by threshold" << LogIO::POST;
2525 232 : timer.mark();
2526 : // make temp mask image consist of the original pix value and below the threshold is set to 0
2527 : //TempImage<Float> maskedRes(res.shape(), res.coordinates(), memoryToUse());
2528 : // single plane
2529 464 : TempImage<Float> maskedRes(planeshp, planeResImage.coordinates(), memoryToUse());
2530 232 : maskedRes.set(0);
2531 : //makeMaskByPerChanThreshold(res, chanFlag, maskedRes, maskThreshold, dummysizes);
2532 232 : makeMaskByPerChanThreshold(planeResImage, chanFlag1elem, maskedRes, maskThreshold1elem, dummysizes);
2533 232 : timeInitThresh(0)+=timer.real(); timeInitThresh(1)+=timer.user(); timeInitThresh(2)+=timer.system();
2534 232 : if (perplanetiming) {
2535 : os << LogIO::NORMAL3 << "End thresholding: time to create the initial threshold mask: real "<< timer.real()
2536 232 : << "s ( user " << timer.user() <<"s, system "<< timer.system() << "s)" << LogIO::POST;
2537 : }
2538 : //if (res.hasPixelMask()) {
2539 232 : if (planeResImage.hasPixelMask()) {
2540 : //os << LogIO::DEBUG1 <<" HAS MASK....ich="<<ich<<LogIO::POST;
2541 232 : ArrayLattice<Bool> pixmasklat(planeResImage.getMask());
2542 232 : maskedRes.copyData( (LatticeExpr<Float>)( iif(pixmasklat, maskedRes, 0.0 ) ) );
2543 : //for debug
2544 : //Array<Float> testdata;
2545 : //maskedRes.get(testdata);
2546 : //os<<" current total of pix values="<<sum(testdata)<<LogIO::POST;
2547 : }
2548 : //TODO MOVE THIS SECTION outside the for-loop -DONE
2549 : //this section need to be move to the end of automask outside of the main chan loop
2550 : //Vector<Bool> allPruned(nchan);
2551 : //if (!iterdone) noMaskCheck(maskedRes, ThresholdType(ipol,ich));
2552 232 : if (debug2) {
2553 0 : os<<LogIO::DEBUG2<<"Saving intermediate masks for this cycle: with name tmp****-"<<iterdone<<".im"<<LogIO::POST;
2554 0 : String tmpfname1="tmpInitThresh-ch"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2555 0 : PagedImage<Float> savedPreMask(planeResImage.shape(),planeResImage.coordinates(),tmpfname1);
2556 0 : savedPreMask.copyData(maskedRes);
2557 : }
2558 :
2559 232 : os << LogIO::NORMAL3 << "Start pruning: the initial threshold mask" << LogIO::POST;
2560 232 : timer.mark();
2561 : //std::shared_ptr<ImageInterface<Float> > tempIm_ptr = pruneRegions2(maskedRes, tempthresh, -1, pruneSize);
2562 : //temporary storage for single plane results
2563 464 : Vector<uInt> nreg1elem, npruned1elem;
2564 464 : Vector<Bool> allPruned1elem(1);
2565 464 : std::shared_ptr<ImageInterface<Float> > tempIm_ptr = YAPruneRegions(maskedRes, chanFlag1elem, allPruned1elem, nreg1elem, npruned1elem, pruneSize, false);
2566 232 : nreg(ipol,ich) = nreg1elem(0);
2567 232 : npruned(ipol, ich) = npruned1elem(0);
2568 232 : allPruned(ipol, ich) = allPruned1elem(0);
2569 : //tempmask.copyData(*(tempIm_ptr.get()));
2570 232 : planeTempMask.copyData(*(tempIm_ptr.get()));
2571 : // now this timing for single plane... need to accumlate and report it later????
2572 232 : timePrune(0)+=timer.real(); timePrune(1)+=timer.user(); timePrune(2)+=timer.system();
2573 232 : if (perplanetiming) {
2574 : os << LogIO::NORMAL3 << "End pruning: time to prune the initial threshold mask: real "
2575 232 : << timer.real()<< "s (user " << timer.user() <<"s, system "<< timer.system() << "s)" << LogIO::POST;
2576 : }
2577 :
2578 : //if (debug2) {
2579 : // String tmpfname2="tmpAfterPrune-"+String::toString(iterdone)+".im";
2580 : // PagedImage<Float> savedPrunedPreThreshMask(res.shape(),res.coordinates(),tmpfname2);
2581 : // savedPrunedPreThreshMask.copyData(*(tempIm_ptr.get()));
2582 : //}
2583 : //themask = LatticeExpr<Float> ( iif( *(tempIm_ptr.get()) > maskThreshold, 1.0, 0.0 ));
2584 : // Need this?
2585 : //makeMaskByPerChanThreshold(*(tempIm_ptr.get()), tempmask, maskThreshold, dummysizes);
2586 : //if (debug) {
2587 : // PagedImage<Float> savedPostPrunedMask(res.shape(),res.coordinates(),"tmp-postPruningPostThreshMask.im");
2588 : // savedPostPrunedMask.copyData(tempmask);
2589 : //}
2590 : }
2591 : else { // ***** No pruning case ******
2592 0 : os << LogIO::NORMAL3 << "Start thresholding: create an initial threshold mask" << LogIO::POST;
2593 0 : timer.mark();
2594 : //tempmask.set(0);
2595 0 : planeTempMask.set(0);
2596 : //makeMaskByPerChanThreshold(res, chanFlag, tempmask, maskThreshold, dummysizes);
2597 0 : makeMaskByPerChanThreshold(planeResImage, chanFlag1elem, planeTempMask, maskThreshold1elem, dummysizes);
2598 0 : if (debug2) {
2599 0 : String tmpnopruneinitmask="tmpInitThreshNoprune-ch"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2600 0 : PagedImage<Float> savedThreshmask(planeResImage.shape(), planeResImage.coordinates(), tmpnopruneinitmask);
2601 0 : savedThreshmask.copyData(planeTempMask);
2602 : }
2603 :
2604 : //if (!iterdone) noMaskCheck(tempmask, ThresholdType);
2605 0 : timePrune(0)+=timer.real(); timePrune(1)+=timer.user(); timePrune(2)+=timer.system();
2606 0 : if (perplanetiming) {
2607 : os << LogIO::NORMAL3 << "End trehsholding: time to create the initial threshold mask: real "
2608 0 : << timer.real()<<"s (user " << timer.user() <<"s, system "<< timer.system() << "s)" << LogIO::POST;
2609 : }
2610 : //tempmask.copyData(themask);
2611 : } // DONE PRUNING STAGE
2612 :
2613 : // ***** SMOOTHING *******
2614 232 : os << LogIO::NORMAL3 << "Start smoothing: the initial threshold mask" << LogIO::POST;
2615 232 : timer.mark();
2616 464 : SPIIF outmask = convolveMask(planeTempMask, modbeam );
2617 232 : if (debug2 ) {
2618 0 : String tmpfname3="tmp-postSmoothMask-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2619 0 : PagedImage<Float> savedSmoothedMask(planeResImage.shape(),planeResImage.coordinates(),tmpfname3);
2620 0 : savedSmoothedMask.copyData(*(outmask.get()));
2621 : }
2622 :
2623 :
2624 : //clean up (appy cutThreshold to convolved mask image)
2625 464 : String lelmask("");
2626 : //use standard stats
2627 464 : Record smmaskstats = calcImageStatistics(*outmask, lelmask, 0, false);
2628 464 : Array<Double> smmaskmaxs;
2629 232 : smmaskstats.get(RecordFieldId("max"),smmaskmaxs);
2630 464 : Vector<Float> cutThresholdValue(1);
2631 :
2632 : //if (npol<=1) {
2633 : // chindx(0) = 0;
2634 : //}
2635 : //else {
2636 : // chindx(1) = 0;
2637 : //}
2638 : // cutThresholdValue(ich) = cutThreshold * smmaskmaxs(chindx);
2639 232 : cutThresholdValue(0) = cutThreshold * smmaskmaxs(IPosition(1,0));
2640 : //os<<LogIO::DEBUG1<<" cutThreshVal("<<ich<<")="<<cutThresholdValue(ich)<<LogIO::POST;
2641 :
2642 : //TempImage<Float> thenewmask(res.shape(),res.coordinates(), memoryToUse());
2643 : //thenewmask.set(0);
2644 464 : TempImage<Float> thenewmask(planeshp,planeResImage.coordinates(), memoryToUse());
2645 232 : thenewmask.set(0);
2646 : //makeMaskByPerChanThreshold(*outmask, chanFlag, thenewmask, cutThresholdValue, dummysizes);
2647 232 : makeMaskByPerChanThreshold(*outmask, chanFlag1elem, thenewmask, cutThresholdValue, dummysizes);
2648 : // Per-plane timing
2649 232 : timeSmooth(0) += timer.real(); timeSmooth(1) += timer.user(); timeSmooth(2) += timer.system();
2650 232 : if (perplanetiming) {
2651 : os << LogIO::NORMAL3 << "End smoothing: time to create the smoothed initial threshold mask: real "<< timer.real()
2652 232 : <<"s (user " << timer.user() <<"s, system "<< timer.system() << "s)" << LogIO::POST;
2653 : }
2654 :
2655 : /***
2656 : if (!iterdone) {
2657 : if (!isEmptyMask(*(outmask.get())) && isEmptyMask(thenewmask)) os<<LogIO::WARN<<"Removed all regions based by cutthreshold applied to the smoothed mask."<<LogIO::POST;
2658 : }
2659 : ***/
2660 232 : if (debug2 ) {
2661 0 : String tmpnewmask="tmp-AfterSmooth-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2662 0 : PagedImage<Float> savedthenewmask(planeResImage.shape(), planeResImage.coordinates(), tmpnewmask);
2663 0 : savedthenewmask.copyData(thenewmask);
2664 : }
2665 :
2666 : // ***** GROW STAGE *****
2667 : //
2668 : // take stats on the current mask for setting flags for grow mask : if max < 1 for any spectral plane it will grow the previous mask
2669 : //
2670 : // Mod: 2017.07.26: modified get stats for prev mask, if channel contains no mask in prev mask it will set flag to skip the channel
2671 : //Record maskstats = calcImageStatistics(thenewmask, thenewmask, lelmask, 0, false);
2672 : //
2673 464 : SubImage<Float> subMask(mask,sl, true, aspec, true);
2674 232 : if(debug2) {
2675 0 : String tmpsubmaskName = "tmp-submask-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2676 0 : PagedImage<Float> tmpsubmask(planeResImage.shape(), planeResImage.coordinates(), tmpsubmaskName);
2677 0 : tmpsubmask.copyData(subMask);
2678 : }
2679 : //Record maskstats = calcImageStatistics(mask, lelmask, 0, false);
2680 : // per-plane stats now
2681 464 : Record maskstats = calcImageStatistics(subMask, lelmask, 0, false);
2682 464 : Array<Double> maskmaxs;
2683 232 : maskstats.get(RecordFieldId("max"),maskmaxs);
2684 464 : IPosition arrshape = maskmaxs.shape();
2685 232 : uInt naxis=arrshape.size();
2686 464 : IPosition indx(naxis,0);
2687 : //os<<LogIO::NORMAL<<"arrshape="<<arrshape<<" indx="<<indx<<LogIO::POST;
2688 : //os<<LogIO::NORMAL<<"statshp="<<statshp<<LogIO::POST;
2689 : // ignoring corr for now and assume first axis is channel
2690 464 : Array<Bool> dogrow(arrshape);
2691 232 : dogrow.set(false);
2692 464 : for (uInt i=0; i < arrshape(0); i++) {
2693 232 : indx(0) = i;
2694 232 : if (maskmaxs(indx) == 1.0 && !chanFlag1elem(0)) {
2695 79 : dogrow(indx) = true;
2696 : }
2697 : //For debug
2698 : //if (chanFlag(i)) {
2699 : // os<<LogIO::NORMAL<<"For dogrow: skipping channel: "<<i<<" chanFlag(i)="<<chanFlag(i)<<" dogrow("<< indx << ")=" <<dogrow(indx)<<LogIO::POST;
2700 : //}
2701 : // set dogrow true for all chans (contraintMask should be able to handle skipping channels )
2702 : // dogrow(indx) = true;
2703 : }
2704 :
2705 232 : if (iterdone && growIterations>0) { // enter to acutal grow process
2706 : //per-plane timing
2707 111 : os << LogIO::NORMAL3 << "Start grow mask: growing the previous mask " << LogIO::POST;
2708 111 : timer.mark();
2709 : //call growMask
2710 : // corresponds to calcThresholdMask with lowNoiseThreshold...
2711 : //TempImage<Float> constraintMaskImage(res.shape(), res.coordinates(), memoryToUse());
2712 : // per-plane constraint mask image
2713 222 : TempImage<Float> constraintMaskImage(planeshp, planeResImage.coordinates(), memoryToUse());
2714 111 : constraintMaskImage.set(0);
2715 : // constrainMask is 1/0 mask
2716 : //makeMaskByPerChanThreshold(res, chanFlag, constraintMaskImage, lowMaskThreshold, dummysizes);
2717 111 : makeMaskByPerChanThreshold(planeResImage, chanFlag1elem, constraintMaskImage, lowMaskThreshold1elem, dummysizes);
2718 111 : if(debug2 && ipol==0 && ich==0) {
2719 0 : os<< LogIO::NORMAL3<<"saving constraint mask " << LogIO::POST;
2720 0 : PagedImage<Float> beforepruneconstIm(planeResImage.shape(), planeResImage.coordinates(),"tmpConstraint-"+String::toString(iterdone)+".im");
2721 0 : beforepruneconstIm.copyData(constraintMaskImage);
2722 : }
2723 :
2724 : // 2017.05.05: should done after multiply by binary dilation
2725 : //
2726 : // prune the constraintImage
2727 : //if (minBeamFrac > 0.0 ) {
2728 : // //Double thethresh=0.1;
2729 : // os<<LogIO::NORMAL << "Pruning the constraint mask "<<LogIO::POST;
2730 : // //std::shared_ptr<ImageInterface<Float> > tempPrunedMask_ptr = pruneRegions2(constraintMaskImage, thethresh, -1, pruneSize);
2731 : // Vector<Bool> dummy(0);
2732 : // std::shared_ptr<ImageInterface<Float> > tempPrunedMask_ptr = YAPruneRegions(constraintMaskImage, dummy, pruneSize);
2733 : // constraintMaskImage.copyData( *(tempPrunedMask_ptr.get()) );
2734 : //}
2735 : //if(debug2) {
2736 : // PagedImage<Float> afterpruneconstIm(res.shape(), res.coordinates(),"tmpAfterPruneConstraint-"+String::toString(iterdone)+".im");
2737 : // afterpruneconstIm.copyData(constraintMaskImage);
2738 : //}
2739 :
2740 : // for mask in binaryDilation, translate it to T/F (if T it will grow the mask region (NOTE currently binary dilation
2741 : // does opposite T/F interpretation NEED to CHANGE)
2742 : //TempImage<Bool> constraintMask(res.shape(),res.coordinates(), memoryToUse());
2743 : //constraintMask.copyData( LatticeExpr<Bool> (iif(constraintMaskImage > 0, true, false)) );
2744 222 : TempImage<Bool> constraintMask(planeshp, planeResImage.coordinates(), memoryToUse());
2745 111 : constraintMask.copyData( LatticeExpr<Bool> (iif(constraintMaskImage > 0, true, false)) );
2746 : // simple structure element for binary dilation
2747 222 : IPosition axislen(2, 3, 3);
2748 222 : Array<Float> se(axislen);
2749 111 : se.set(0);
2750 111 : se(IPosition(2,1,0))=1.0;
2751 111 : se(IPosition(2,0,1))=1.0;
2752 111 : se(IPosition(2,1,1))=1.0;
2753 111 : se(IPosition(2,2,1))=1.0;
2754 111 : se(IPosition(2,1,2))=1.0;
2755 111 : if(debug2 && ich==0 && ipol == 0) {
2756 0 : String tmpbeforeBD="tmp-BeforeBD-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2757 0 : PagedImage<Float> beforeBinaryDilationIm(planeResImage.shape(), planeResImage.coordinates(),tmpbeforeBD);
2758 : //`beforeBinaryDilationIm.copyData(constraintMaskImage);
2759 0 : beforeBinaryDilationIm.copyData(subMask);
2760 : }
2761 : // CHECK THIS works for a single plane but 4 dim image
2762 :
2763 111 : subprevmask.set(0);
2764 111 : binaryDilation(subMask, se, growIterations, constraintMask, dogrow, subprevmask);
2765 111 : if(debug2) {
2766 0 : PagedImage<Float> afterBinaryDilationIm(planeResImage.shape(), planeResImage.coordinates(),"tmpAfterBinaryDilation-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im");
2767 0 : afterBinaryDilationIm.copyData(subprevmask);
2768 : }
2769 : // multiply binary dilated mask by constraintmask
2770 : //prevmask.copyData( LatticeExpr<Float> (constraintMaskImage*prevmask));
2771 111 : subprevmask.copyData( LatticeExpr<Float> (constraintMaskImage*subprevmask));
2772 111 : if(debug2) {
2773 0 : PagedImage<Float> beforepruneconstIm(planeResImage.shape(), planeResImage.coordinates(),"tmpBeforePruneGrowMask-"+String::toString(ich)+"ipol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im");
2774 0 : beforepruneconstIm.copyData(subprevmask);
2775 : }
2776 111 : timeGrow(0) += timer.real(); timeGrow(1) += timer.user(); timeGrow(2) += timer.system();
2777 111 : if (perplanetiming) {
2778 : os << LogIO::NORMAL3 << "End grow mask: time to grow the previous mask: real "
2779 111 : << timer.real() <<"s (user "<< timer.user() << "s, system " << timer.system() << "s)" << LogIO::POST;
2780 : }
2781 :
2782 : // **** pruning on grow mask ****
2783 111 : if (minBeamFrac > 0.0 && doGrowPrune) {
2784 : //os<<LogIO::NORMAL << "Pruning the growed previous mask "<<LogIO::POST;
2785 111 : os << LogIO::NORMAL3 << "Start pruning: on the grow mask "<< LogIO::POST;
2786 111 : timer.mark();
2787 222 : Vector<Bool> dummy(0);
2788 222 : Vector<uInt> ngrowreg1elem, ngrowpruned1elem;
2789 : //std::shared_ptr<ImageInterface<Float> > tempPrunedMask_ptr = YAPruneRegions(prevmask, chanFlag, dummy, ngrowreg, ngrowpruned, pruneSize);
2790 :
2791 222 : std::shared_ptr<ImageInterface<Float> > tempPrunedMask_ptr = YAPruneRegions(subprevmask, chanFlag1elem, dummy, ngrowreg1elem, ngrowpruned1elem, pruneSize, false);
2792 111 : ngrowreg(ipol,ich) = ngrowreg1elem(0);
2793 111 : ngrowpruned(ipol,ich) = ngrowpruned1elem(0);
2794 : //prevmask.copyData( *(tempPrunedMask_ptr.get()) );
2795 111 : subprevmask.copyData( *(tempPrunedMask_ptr.get()) );
2796 111 : timePruneGrow(0) += timer.real(); timePruneGrow(1) += timer.user(); timePruneGrow(2) += timer.system();
2797 111 : if (perplanetiming) {
2798 : os << LogIO::NORMAL3 << "End pruning: time to prune the grow mask: real "
2799 111 : << timer.real() <<"s (user "<< timer.user() << "s, system "<< timer.system() << "s)" << LogIO::POST;
2800 : }
2801 111 : if(debug2) {
2802 0 : String tmpafterprunegrowname="tmpAfterPruneGrowMask-"+String::toString(ich)+"pol"+String::toString(ipol)+"iter"+String::toString(iterdone)+".im";
2803 0 : PagedImage<Float> afterpruneconstIm(planeResImage.shape(), planeResImage.coordinates(),tmpafterprunegrowname);
2804 0 : afterpruneconstIm.copyData(subprevmask);
2805 : }
2806 : }
2807 :
2808 : // ***** smoothing on grow mask *****
2809 111 : os << LogIO::NORMAL3 << "Start smoothing: the grow mask " << LogIO::POST;
2810 111 : timer.mark();
2811 : ///SPIIF outprevmask = convolveMask( prevmask, modbeam);
2812 222 : SPIIF outprevmask = convolveMask( subprevmask, modbeam);
2813 : //if (debug) {
2814 : // PagedImage<Float> postSmoothGrowedMask(res.shape(), res.coordinates(),"tmpPostSmoothGrowMask-"+String::toString(iterdone)+".im");
2815 : // postSmoothGrowedMask.copyData(*outprevmask);
2816 : //}
2817 : //prevmask.copyData( LatticeExpr<Float> (iif( *(outprevmask.get()) > cutThreshold, 1.0, 0.0 )) );
2818 : // single plane
2819 222 : Record constmaskstats = calcImageStatistics(*outprevmask, lelmask, 0, false);
2820 222 : Array<Double> constmaskmaxs;
2821 111 : constmaskstats.get(RecordFieldId("max"),constmaskmaxs);
2822 : //Vector<Float> constCutThresholdValue(nchan);
2823 222 : Vector<Float> constCutThresholdValue(1);
2824 : // stats on a single plane now, so no need of chindx
2825 : //if (npol <=1) {
2826 : // chindx(0) = 0;
2827 : //}
2828 : //else {
2829 : // chindx(1) = 0;
2830 : //}
2831 : //constCutThresholdValue(0) = cutThreshold * constmaskmaxs(chindx);
2832 111 : constCutThresholdValue(0) = cutThreshold * constmaskmaxs(IPosition(1,0));
2833 : //prevmask.set(0);
2834 111 : subprevmask.set(0);
2835 : //makeMaskByPerChanThreshold(*outprevmask, chanFlag, prevmask, constCutThresholdValue, dummysizes);
2836 111 : makeMaskByPerChanThreshold(*outprevmask, chanFlag1elem, subprevmask, constCutThresholdValue, dummysizes);
2837 : //if (debug) {
2838 : // PagedImage<Float> smoothedGrowedMask(res.shape(), res.coordinates(),"tmpSmoothedGrowMask-"+String::toString(iterdone)+".im");
2839 : // smoothedGrowedMask.copyData(prevmask);
2840 : //}
2841 111 : timeSmoothGrow(0) += timer.real(); timeSmoothGrow(1) += timer.user(); timeSmoothGrow(2) += timer.system();
2842 111 : if (perplanetiming) {
2843 : os << LogIO::NORMAL3 << "End smoothing: time to create the smoothed grow mask: real "
2844 111 : << timer.real() <<"s (user "<< timer.user() << "s, system " << timer.system() << "s)" << LogIO::POST;
2845 : }
2846 : } //end - GROW (iterdone && dogrowiteration)
2847 :
2848 : // ****** save positive (emission) mask only ******
2849 :
2850 : // temporary save negative mask from the previous one
2851 : //TempImage<Float> prevnegmask(res.shape(), res.coordinates(), memoryToUse());
2852 : //prevnegmask.copyData( (LatticeExpr<Float>)( iif( (mask - posmask ) > 0.0, 1.0, 0.0 ) ) );
2853 :
2854 : //
2855 232 : if (debug2 ) {
2856 0 : String beforesumSPmaskname= "beforesumSPmask-ch"+String::toString(ich)+"pol"+String::toString(ipol)+"-iter"+String::toString(iterdone)+".im";
2857 0 : PagedImage<Float> tempsubposmask(TiledShape(subposmask.shape()), subposmask.coordinates(), beforesumSPmaskname);
2858 0 : tempsubposmask.copyData(subposmask);
2859 0 : String beforesumSPrevmaskname = "beforesumSPrevmask-ch"+String::toString(ich)+"pol"+String::toString(ipol)+"-iter"+String::toString(iterdone)+".im";
2860 0 : PagedImage<Float> tempsubprevmask(TiledShape(subposmask.shape()), subposmask.coordinates(), beforesumSPrevmaskname);
2861 0 : tempsubprevmask.copyData(subprevmask);
2862 : }
2863 :
2864 :
2865 :
2866 232 : if (planeResImage.hasPixelMask()) {
2867 : //LatticeExpr<Bool> pixmask(res.pixelMask());
2868 464 : LatticeExpr<Bool> pixmask(planeResImage.pixelMask());
2869 : // add all positive masks (previous one, grow mask, current thresh mask)
2870 : // mask = untouched prev mask, prevmask=modified prev mask by the grow func, thenewmask=mask by thresh on current residual
2871 : //posmask.copyData( (LatticeExpr<Float>)( iif((posmask + prevmask + thenewmask ) > 0.0 && pixmask, 1.0, 0.0 ) ) );
2872 232 : subposmask.copyData( (LatticeExpr<Float>)( iif((subposmask + subprevmask + thenewmask ) > 0.0 && pixmask, 1.0, 0.0 ) ) );
2873 232 : os <<LogIO::DEBUG1 <<"Add positive previous mask, pbmask and the new mask for this plane"<<LogIO::POST;
2874 : }
2875 : else {
2876 : //posmask.copyData( (LatticeExpr<Float>)( iif((posmask + prevmask + thenewmask ) > 0.0, 1.0, 0.0 ) ) );
2877 0 : subposmask.copyData( (LatticeExpr<Float>)( iif((subposmask + subprevmask + thenewmask ) > 0.0, 1.0, 0.0 ) ) );
2878 0 : os <<LogIO::DEBUG1 <<"Add positive previous mask and the new mask.."<<LogIO::POST;
2879 : }
2880 :
2881 : // **** NEGATIVE MASK creation *****
2882 : //TempImage<Float> thenegmask(res.shape(),res.coordinates(), memoryToUse());
2883 464 : TempImage<Float> subnegmask(planeshp, planeResImage.coordinates(), memoryToUse());
2884 232 : subnegmask.set(0);
2885 : //Vector<Float> negmaskpixs;
2886 464 : Vector<Float> negmaskpixs1elem;
2887 232 : if (negativeThresholdFactor > 0) {
2888 5 : os << LogIO::NORMAL3 << "Start thresholding: create a negative mask" << LogIO::POST;
2889 5 : timer.mark();
2890 : //os<<LogIO::NORMAL<<"Creating a mask for negative features. "<<LogIO::POST;
2891 : //TempImage<Float> negativeMaskImage(res.shape(), res.coordinates(), memoryToUse());
2892 10 : TempImage<Float> negativeSubMaskImage(planeshp, planeResImage.coordinates(), memoryToUse());
2893 5 : negativeSubMaskImage.set(0);
2894 : //makeMaskByPerChanThreshold(res, chanFlag, negativeMaskImage , negativeMaskThreshold, dummysizes);
2895 5 : makeMaskByPerChanThreshold(planeResImage, chanFlag1elem, negativeSubMaskImage, negativeMaskThreshold1elem, dummysizes);
2896 : // SPIIF negmask = convolveMask( negativeMaskImage, modbeam);
2897 10 : SPIIF negmask = convolveMask( negativeSubMaskImage, modbeam);
2898 : // determine the cutthreshold value for negative mask
2899 10 : Record negmaskstats = calcImageStatistics(*negmask, lelmask, 0, false);
2900 10 : Array<Double> negmaskmaxs;
2901 5 : negmaskstats.get(RecordFieldId("max"),negmaskmaxs);
2902 : //Vector<Float> negCutThresholdValue(nchan);
2903 10 : Vector<Float> negCutThresholdValue(1);
2904 : // 1 dim stats now, so no need of chindx
2905 : //if (npol <= 1) {
2906 : // chindx(0) = 0;
2907 : //}
2908 : //else {
2909 : // chindx(1) = 0;
2910 : //}
2911 5 : negCutThresholdValue(0) = cutThreshold * negmaskmaxs(IPosition(1,0));
2912 : //makeMaskByPerChanThreshold(*negmask, chanFlag, thenegmask, negCutThresholdValue, negmaskpixs);
2913 5 : makeMaskByPerChanThreshold(*negmask, chanFlag1elem, subnegmask, negCutThresholdValue, negmaskpixs1elem);
2914 5 : negmaskpixs(ipol,ich) = negmaskpixs1elem(0);
2915 5 : if (isEmptyMask(subnegmask) ){
2916 1 : os<<"No negative region was found by auotmask."<<LogIO::POST;
2917 : }
2918 : //if (debug) {
2919 : // PagedImage<Float> temppresmoothnegmask(TiledShape(negativeMaskImage.shape()), negativeMaskImage.coordinates(),"tmpPreSmoNegMask.im");
2920 : // temppresmoothnegmask.copyData(negativeMaskImage);
2921 : //PagedImage<Float> tempnegmask(TiledShape(thenegmask.shape()), thenegmask.coordinates(),"tmpNegMask.im");
2922 : // tempnegmask.copyData(thenegmask);
2923 : // PagedImage<Float> tempsmonegmask(TiledShape(thenegmask.shape()), thenegmask.coordinates(),"tmpSmoNegMask.im");
2924 : // tempsmonegmask.copyData(*negmask);
2925 : //}
2926 5 : timeNegThresh(0) += timer.real(); timeNegThresh(1) += timer.user(); timeNegThresh(2) += timer.system();
2927 5 : if (perplanetiming) {
2928 : os << LogIO::NORMAL3 << "End thresholding: time to create the negative mask: real "
2929 5 : << timer.real() <<"s (user " << timer.user() << "s, system " << timer.system() << "s)" << LogIO::POST;
2930 : }
2931 : }
2932 :
2933 : // store per plane masks to full cube mask images
2934 : // subMask to mask, posmask, thenegmask
2935 464 : Array<Float> maskdata, posmaskdata, negmaskdata;
2936 : // all images are single planes
2937 464 : IPosition stride(4,1,1,1,1);
2938 464 : IPosition plstart(4,0,0,0,0);
2939 464 : Slicer plsl(plstart, length);
2940 232 : subMask.doGetSlice(maskdata, plsl);
2941 232 : subposmask.doGetSlice(posmaskdata, plsl);
2942 232 : subnegmask.getSlice(negmaskdata, plsl);
2943 : // put to full mask images
2944 232 : mask.putSlice(maskdata,start,stride);
2945 232 : posmask.putSlice(posmaskdata,start,stride);
2946 232 : thenegmask.putSlice(negmaskdata,start,stride);
2947 232 : if (debug2) {
2948 0 : PagedImage<Float> tempfinalmask(TiledShape(mask.shape()), mask.coordinates(), "tmpfinalmask"+String::toString(iterdone)+".im");
2949 0 : tempfinalmask.copyData(mask);
2950 : }
2951 : } // if not chanFlag=True
2952 : } // ipol iter
2953 : } // the main per plane for-loop end for-ich
2954 :
2955 : //print tot. timing for each step
2956 52 : if (npol*nchan > 1) {
2957 12 : os << LogIO::NORMAL << "*** Timing summary for whole planes ***" << LogIO::POST;
2958 24 : os << LogIO::NORMAL << "Total time to create the initial threshold mask: real "<< timeInitThresh(0)
2959 12 : << "s ( user " << timeInitThresh(1) <<"s, system "<< timeInitThresh(2) << "s)" << LogIO::POST;
2960 : os << LogIO::NORMAL << "Total time to prune the initial threshold mask: real "
2961 12 : << timePrune(0)<< "s (user " << timePrune(1) <<"s, system "<< timePrune(2) << "s)" << LogIO::POST;
2962 24 : os << LogIO::NORMAL << "Total time to create the smoothed initial threshold mask: real "<< timeSmooth(0)
2963 12 : <<"s (user " << timeSmooth(1)<<"s, system "<< timeSmooth(2) << "s)" << LogIO::POST;
2964 : os << LogIO::NORMAL << "Total time to grow the previous mask: real "
2965 12 : << timeGrow(0) <<"s (user "<< timeGrow(1) << "s, system " << timeGrow(2) << "s)" << LogIO::POST;
2966 : os << LogIO::NORMAL << "Total time to prune the grow mask: real "
2967 12 : << timePruneGrow(0) <<"s (user "<< timePruneGrow(1) << "s, system "<< timePruneGrow(2) << "s)" << LogIO::POST;
2968 : os << LogIO::NORMAL << "Total time to create the smoothed grow mask: real "
2969 12 : << timeSmoothGrow(0) <<"s (user "<< timeSmoothGrow(1) << "s, system " << timeSmoothGrow(2) << "s)" << LogIO::POST;
2970 : }
2971 :
2972 52 : if (!iterdone) noMaskCheck(posmask, ThresholdType);
2973 : // "Allpruned check here"
2974 52 : Int nAllPruned=ntrue(allPruned);
2975 : // pruning is done only on positive mask
2976 52 : if(!iterdone && isEmptyMask(posmask) && nAllPruned) {
2977 : os<<LogIO::WARN<<nAllPruned<<" of "<<nchan<<" channels had all regions removed by pruning."
2978 0 : <<" Try decreasing minbeamfrac to remove fewer regions"<<LogIO::POST;
2979 : }
2980 :
2981 : //for debug
2982 : /***
2983 : PagedImage<Float> tempthemask(TiledShape(tempIm_ptr.get()->shape()), tempIm_ptr.get()->coordinates(),"tempthemask.Im");
2984 : tempthemask.copyData(themask);
2985 : ***/
2986 :
2987 : // In the initial iteration, if no mask is created (all spectral planes) by automask it will fall back to full clean mask
2988 : /***
2989 : if (!iterdone) {
2990 : Array<Float> maskdata;
2991 : IPosition maskshape = thenewmask.shape();
2992 : Int naxis = maskshape.size();
2993 : IPosition blc(naxis,0);
2994 : IPosition trc=maskshape-1;
2995 : Slicer sl(blc,trc,Slicer::endIsLast);
2996 : thenewmask.doGetSlice(maskdata,sl);
2997 : if (sum(maskdata)==0.0) {
2998 : mask.set(1);
2999 : //os<<LogIO::WARN<<"No mask was created by automask, set a clean mask to the entire image."<<LogIO::POST;
3000 : os<<LogIO::WARN<<"No mask was created by automask."<<LogIO::POST;
3001 : }
3002 : }
3003 : ***/
3004 52 : if (debug) {
3005 : //saved prev unmodified mask
3006 0 : PagedImage<Float> tmpUntouchedPrevMask(res.shape(), res.coordinates(),"tmpUnmodPrevMask"+String::toString(iterdone)+".im");
3007 0 : tmpUntouchedPrevMask.copyData(mask);
3008 :
3009 : }
3010 : // make a copy of unmodified previous mask
3011 156 : TempImage<Float> unmodifiedprevmask(res.shape(),res.coordinates(), memoryToUse());
3012 52 : unmodifiedprevmask.copyData(mask);
3013 :
3014 52 : if (res.hasPixelMask()) {
3015 104 : LatticeExpr<Bool> pixmask(res.pixelMask());
3016 : //mask.copyData( (LatticeExpr<Float>)( iif((mask + thenewmask) > 0.0 && pixmask, 1.0, 0.0 ) ) );
3017 : // add all masks (previous one, grow mask, current thresh mask)
3018 : // mask = untouched prev mask, prevmask=modified prev mask by the grow func, thenewmask=mask by thresh on current residual
3019 : //mask.copyData( (LatticeExpr<Float>)( iif((mask+prevmask + thenewmask + thenegmask) > 0.0 && pixmask, 1.0, 0.0 ) ) );
3020 52 : if(debug) {
3021 0 : PagedImage<Float> savedUnmod(res.shape(), res.coordinates(), "savedUmmod"+String::toString(iterdone)+".im");
3022 0 : savedUnmod.copyData(mask);
3023 0 : PagedImage<Float> savedPosmask(res.shape(), res.coordinates(), "savedPosmask"+String::toString(iterdone)+".im");
3024 0 : savedPosmask.copyData(posmask);
3025 0 : PagedImage<Float> savedNegmask(res.shape(), res.coordinates(), "savedNegmask"+String::toString(iterdone)+".im");
3026 0 : savedPosmask.copyData(thenegmask);
3027 : }
3028 52 : mask.copyData( (LatticeExpr<Float>)( iif((mask + posmask + thenegmask ) > 0.0 && pixmask, 1.0, 0.0 ) ) );
3029 :
3030 52 : mask.clearCache();
3031 52 : mask.unlock();
3032 52 : mask.tempClose();
3033 52 : os <<LogIO::DEBUG1 <<"Add previous mask, pbmask and the new mask.."<<LogIO::POST;
3034 : }
3035 : else {
3036 : //os <<"Lattice themask is created..."<<LogIO::POST;
3037 : //LatticeExpr<Float> themask( iif( tempconvim > rmsthresh/afactor, 1.0, 0.0 ));
3038 :
3039 : //mask.copyData( (LatticeExpr<Float>)( iif((mask + prevmask + thenewmask + thenegmask ) > 0.0, 1.0, 0.0 ) ) );
3040 0 : mask.copyData( (LatticeExpr<Float>)( iif((mask + posmask + thenegmask ) > 0.0, 1.0, 0.0 ) ) );
3041 :
3042 0 : os <<LogIO::DEBUG1 <<"Add previous mask and the new mask.."<<LogIO::POST;
3043 : }
3044 : // test the curent final mask with the previous mask
3045 104 : Vector<Bool> zeroChanMask;
3046 52 : skipChannels(fracChange,unmodifiedprevmask, mask, ThresholdType, isthresholdreached, chanFlag, zeroChanMask);
3047 :
3048 52 : if (verbose)
3049 4 : printAutomaskSummary(resRmss, maxs, mins, mdns, maskThreshold, ThresholdType, chanFlag, zeroChanMask, nreg, npruned, ngrowreg, ngrowpruned, negmaskpixs, summaryRec);
3050 52 : }//end of autoMaskByMultiThreshold
3051 :
3052 186 : Bool SDMaskHandler::isEmptyMask(ImageInterface<Float>& mask)
3053 : {
3054 372 : Array<Float> maskdata;
3055 372 : IPosition maskshape = mask.shape();
3056 186 : Int naxis = maskshape.size();
3057 372 : IPosition blc(naxis,0);
3058 372 : IPosition trc=maskshape-1;
3059 186 : Slicer sl(blc,trc,Slicer::endIsLast);
3060 186 : mask.doGetSlice(maskdata,sl);
3061 372 : return (sum(maskdata)==0.0);
3062 :
3063 : }
3064 :
3065 : //void SDMaskHandler::noMaskCheck(ImageInterface<Float>& mask, Vector<String>& thresholdType)
3066 31 : void SDMaskHandler::noMaskCheck(ImageInterface<Float>& mask, Matrix<String>& thresholdType)
3067 : {
3068 : // checkType and thresholdType will determine the exact messages to print out in the log
3069 93 : LogIO os( LogOrigin("SDMaskHandler","autoMaskByMultiThreshold",WHERE) );
3070 : // for waring messsages
3071 : /***
3072 : Array<Float> maskdata;
3073 : IPosition maskshape = mask.shape();
3074 : Int naxis = maskshape.size();
3075 : IPosition blc(naxis,0);
3076 : IPosition trc=maskshape-1;
3077 : Slicer sl(blc,trc,Slicer::endIsLast);
3078 : mask.doGetSlice(maskdata,sl);
3079 : ***/
3080 : //if (sum(maskdata)==0.0) {
3081 31 : if (isEmptyMask(mask)) {
3082 1 : os << LogIO::WARN <<"No regions found by automasking" <<LogIO::POST;
3083 : // checktype
3084 : //Int nThresholdType = thresholdType.nelements();
3085 1 : Int nrow = thresholdType.nrow();
3086 1 : Int ncol = thresholdType.ncolumn();
3087 : //Int nThresholdType = thresholdType.nelements();
3088 1 : Int nThresholdType = nrow*ncol;
3089 1 : Int nsidelobethresh=0;
3090 1 : Int nnoisethresh=0;
3091 1 : if (nThresholdType>1 ) {
3092 : //for (uInt j=0; j<(uInt) nThresholdType; j++) {
3093 0 : for (uInt j=0; j<(uInt) ncol; j++) {
3094 0 : for (uInt i=0; i<(uInt) nrow; i++) {
3095 0 : if (thresholdType(i,j)=="sidelobe") {
3096 0 : nsidelobethresh++;
3097 : }
3098 0 : if (thresholdType(i,j)=="noise") {
3099 0 : nnoisethresh++;
3100 : }
3101 : }
3102 : }
3103 0 : if (nsidelobethresh) {
3104 : os << LogIO::WARN <<nsidelobethresh<<" of "<<nThresholdType
3105 : <<" channels used the sidelobe threshold to mask, but no emission was found."
3106 0 : <<" Try decreasing your sidelobethreshold parameter if you want to capture emission in these channels."<< LogIO::POST;
3107 : }
3108 0 : if (nnoisethresh) {
3109 : os << LogIO::WARN <<nnoisethresh<<" of "<<nThresholdType<<" channels used the noise threshold to mask, but no emission was found."
3110 0 : << " Try decreasing your noisethreshold parameter if you want to capture emission in these channels."<< LogIO::POST;
3111 : }
3112 : }
3113 : else {
3114 :
3115 2 : os << LogIO::WARN << "Used "<<thresholdType(0,0)<<" threshold to mask, but no emission was found."
3116 1 : << "Try decreasing your "<<thresholdType(0,0)<<"threshold parameter if you want to capture emission in these channels."<< LogIO::POST;
3117 : }
3118 : }
3119 31 : }
3120 :
3121 52 : void SDMaskHandler::skipChannels(const Float& fracChange,
3122 : ImageInterface<Float>& prevmask,
3123 : ImageInterface<Float>& curmask,
3124 : //const Vector<String>& thresholdtype,
3125 : const Matrix<String>& thresholdtype,
3126 : const Bool isthresholdreached,
3127 : Vector<Bool>& chanFlag,
3128 : Vector<Bool>& zeroChanMask)
3129 : {
3130 156 : LogIO os( LogOrigin("SDMaskHandler","skipChannels",WHERE) );
3131 : // debug
3132 52 : os<<LogIO::DEBUG1<<"Inside skipChannels...."<<LogIO::POST;
3133 104 : IPosition shp = curmask.shape();
3134 52 : Int naxis = shp.size();
3135 104 : CoordinateSystem csys = curmask.coordinates();
3136 52 : Int specaxis = CoordinateUtil::findSpectralAxis(csys);
3137 52 : Int nchan = shp(specaxis);
3138 : // Assumption here is skipChannels applied to stokes I only to keep track which channels
3139 : // to skip for new automasking
3140 104 : IPosition blc(naxis,0);
3141 104 : IPosition trc=shp-1;
3142 52 : zeroChanMask.resize(nchan);
3143 52 : os<<LogIO::DEBUG1<<"Inside skipChannels....after zeroChanMask init"<<LogIO::POST;
3144 260 : for (Int ichan=0; ichan<nchan; ichan++) {
3145 208 : blc(specaxis)=ichan;
3146 208 : trc(specaxis)=ichan;
3147 416 : Slicer sl(blc,trc,Slicer::endIsLast);
3148 416 : Array<Float> curmaskdata;
3149 208 : curmask.doGetSlice(curmaskdata,sl);
3150 208 : Float curmaskpix = sum(curmaskdata);
3151 : // sepearately store zero channel mask info maybe combined in future to streamline
3152 208 : if (curmaskpix==0) {
3153 65 : zeroChanMask(ichan) = True;
3154 : }
3155 : else {
3156 143 : zeroChanMask(ichan) = False;
3157 : }
3158 :
3159 : //if (thresholdtype(ichan).contains("noise") && isthresholdreached && !chanFlag(ichan)) {
3160 208 : if (thresholdtype(0, ichan).contains("noise") && !chanFlag(ichan)) {
3161 8 : Array<Float> prevmaskdata;
3162 4 : prevmask.doGetSlice(prevmaskdata,sl);
3163 4 : Float prevmaskpix = sum(prevmaskdata);
3164 : //cerr<<"prevmaskpix="<<prevmaskpix<<" curemaskpix="<<curmaskpix<<endl;
3165 : //cerr<<"fracChnage="<<fracChange<<endl;
3166 4 : Float diffpix = abs(curmaskpix-prevmaskpix);
3167 : // stopmask is true if one of the followings is satified
3168 : // 1) if current mask is zero (curmaskpix==0.0)
3169 : // 2) if cyclethreshold==threshold (i.e. isthresholdreached=True) and diffpix is zero or
3170 : // less than user-specified fractinal change
3171 : //if ( curmaskpix==0.0 || (diffpix == 0.0 && prevmaskpix!=0.0) || diffpix < fracChange*prevmaskpix) {
3172 : //if ( curmaskpix==0.0 || (isthresholdreached && ((diffpix == 0.0 && prevmaskpix!=0.0) || diffpix < fracChange*prevmaskpix)) ) {
3173 4 : if ( curmaskpix==0.0 ||
3174 4 : (fracChange >=0.0 && isthresholdreached && ( diffpix == 0.0 || diffpix < fracChange*prevmaskpix) ) ) {
3175 0 : chanFlag(ichan) = True;
3176 0 : os<<LogIO::NORMAL<<"Stopping masking for chan="<<ichan<<LogIO::POST;
3177 : }
3178 : }
3179 : } // for loop end
3180 52 : }
3181 :
3182 0 : std::shared_ptr<ImageInterface<Float> > SDMaskHandler::makeMaskFromBinnedImage(const ImageInterface<Float>& image,
3183 : const Int nx,
3184 : const Int ny,
3185 : const Float& fracofpeak,
3186 : const Float& sigma,
3187 : const Int nmask,
3188 : const Bool autoadjust,
3189 : Double thresh)
3190 : {
3191 0 : Bool debug(False);
3192 0 : Bool firstrun(False);
3193 0 : LogIO os( LogOrigin("SDMaskHandler","makeMaskfromBinnedImage",WHERE) );
3194 0 : RebinImage<Float> tempRebinnedIm( image, IPosition(4,nx, ny,1,1) );
3195 : // for debug
3196 0 : if (debug) {
3197 0 : PagedImage<Float> copyRebinnedIm(TiledShape(tempRebinnedIm.shape()), tempRebinnedIm.coordinates(), "binned.Im");
3198 0 : copyRebinnedIm.copyData(tempRebinnedIm);
3199 : }
3200 :
3201 : // modified threshold
3202 : // original algortihm
3203 : //thresh = 3.0*thresh / sqrt(npix);
3204 : // modified by bin size only
3205 : //thresh = thresh / sqrt(nx);
3206 :
3207 : //stats on the binned image (the info not used for mask criteria yet)
3208 0 : Array<Double> rmses, maxes, mins;
3209 : // vars to store min,max values of extrema and rms in all planes
3210 : Double minRmsVal, maxRmsVal, minMaxVal, maxMaxVal, minMinVal, maxMinVal;
3211 0 : IPosition minRmsPos, maxRmsPos, minMaxPos, maxMaxPos, minMinPos, maxMinPos;
3212 0 : TempImage<Float>* tempImForStat = new TempImage<Float>(tempRebinnedIm.shape(), tempRebinnedIm.coordinates(), memoryToUse() );
3213 0 : tempImForStat->copyData(tempRebinnedIm);
3214 0 : std::shared_ptr<casacore::ImageInterface<float> > temprebin_ptr(tempImForStat);
3215 : //os<<" temprebin_ptr.get()->hasPixelMask()="<<temprebin_ptr.get()->hasPixelMask()<<LogIO::POST;
3216 0 : ImageStatsCalculator<Float> imcalc( temprebin_ptr, 0, "", False);
3217 0 : Vector<Int> stataxes(2);
3218 0 : stataxes[0] = 0;
3219 0 : stataxes[1] = 1;
3220 0 : imcalc.setAxes(stataxes);
3221 0 : Record imstats = imcalc.statistics();
3222 0 : imstats.get(RecordFieldId("rms"),rmses);
3223 0 : imstats.get(RecordFieldId("max"),maxes);
3224 0 : imstats.get(RecordFieldId("min"),mins);
3225 0 : minMax(minRmsVal,maxRmsVal,minRmsPos, maxRmsPos,rmses);
3226 0 : minMax(minMaxVal,maxMaxVal,minMaxPos, maxMaxPos,maxes);
3227 0 : minMax(minMinVal,maxMinVal,minMinPos, maxMinPos,mins);
3228 :
3229 :
3230 : //os << LogIO::NORMAL <<"Statistics on binned image: max="<<maxMaxVal<<" rms="<<maxRmsVal<<LogIO::POST;
3231 : os << LogIO::NORMAL <<"Statistics on binned image: Peak (max)="<<maxMaxVal<<"@"<<maxMaxPos
3232 0 : <<"rms (max) ="<<maxRmsVal<<"@"<<maxRmsPos<<LogIO::POST;
3233 : //os << LogIO::NORMAL <<"Statistics on binned image: min of Max="<<minMaxVal<<"@"<<minMaxPos<<LogIO::POST;
3234 : //os << LogIO::NORMAL <<"Statistics on binned image: min of rms="<<minRmsVal<<"@"<<minRmsPos<<LogIO::POST;
3235 :
3236 0 : TempImage<Float>* tempMask = new TempImage<Float> (tempRebinnedIm.shape(), tempRebinnedIm.coordinates(), memoryToUse() );
3237 :
3238 :
3239 : //Double dr = maxMaxVal/rmses(maxMaxPos);
3240 :
3241 : // save only the first time
3242 0 : if (itsMax==DBL_MAX) {
3243 0 : itsMax = maxMaxVal;
3244 0 : firstrun = True;
3245 : }
3246 : else {
3247 0 : firstrun = False;
3248 : }
3249 :
3250 0 : Float adjFact = 0.0;
3251 :
3252 : //Float fracDiffMax = (itsMax - maxMaxVal)/itsMax;
3253 0 : Float fracDiffRms = (itsRms - maxRmsVal)/itsRms;
3254 : //os<<"fractional changes in max:"<<fracDiffMax<<" in rms:"<<fracDiffRms<<LogIO::POST;
3255 0 : os<<LogIO::DEBUG1<<"fractional changes in rms from previous one:"<<fracDiffRms<<LogIO::POST;
3256 0 : if (autoadjust) {
3257 : //os <<"autoAdjust is on "<<LogIO::POST;
3258 : // automatically adjust threshold for the initial round and when fractional rms change is
3259 : // is less than 10%
3260 0 : if (fracDiffRms < 0.1 ) {
3261 0 : adjFact = (Float) Int(log(1.0/abs(fracDiffRms))-1.0);
3262 : }
3263 : // else if (dr < 10.0) {
3264 : // os<<LogIO::DEBUG1<<"dynamic range:max/rms = "<<dr<<LogIO::POST;
3265 : // adjFact = sigma!=0 && sigma <= 3.0? 2: 0;
3266 : // }
3267 0 : else if (firstrun) {
3268 0 : adjFact = 3;
3269 : }
3270 : }
3271 0 : if (fracofpeak) {
3272 0 : thresh = fracofpeak * maxMaxVal;
3273 0 : Double prevthresh = thresh;
3274 0 : if (adjFact >0.0 ) {
3275 0 : thresh = max((adjFact + 3) * maxRmsVal,thresh);
3276 0 : if (firstrun) {
3277 : // adjustment at 1st iteration cycle, if threshold is too big, make cutoff at 50% peak
3278 0 : if (thresh < itsMax) {
3279 0 : if (prevthresh != thresh) {
3280 0 : os << LogIO::NORMAL <<"first iteration automatically adjusted thresh="<<thresh<<"( "<<" * (adj fact.="<<adjFact+3<<") * rms )"<<LogIO::POST;
3281 : }
3282 : }
3283 : else {
3284 0 : thresh=prevthresh;
3285 : }
3286 : }
3287 : else {
3288 0 : if (prevthresh != thresh) {
3289 0 : os << LogIO::NORMAL <<"thresh="<<thresh<<" ( adj fact.="<<adjFact+3<<") * rms )"<<LogIO::POST;
3290 : }
3291 : }
3292 : }
3293 : // if sidelobe level is set and if it is larger than the current thresh use that value
3294 : //thresh = max( itsSidelobeLevel*itsMax, thresh );
3295 0 : if (adjFact != 0.0) {
3296 : }
3297 : else {
3298 0 : os << LogIO::NORMAL <<"thresh="<<thresh<<" ( "<<fracofpeak<<"* max peak )"<<LogIO::POST;
3299 : }
3300 : }
3301 0 : else if (sigma) {
3302 0 : thresh = (sigma + adjFact)* maxRmsVal;
3303 : // if sidelobe level is set and if it is larger than the current thresh use that value
3304 : //thresh = max( itsSidelobeLevel*itsMax, thresh);
3305 0 : if (firstrun && adjFact != 0.0) {
3306 0 : if (thresh < itsMax) {
3307 : os << LogIO::NORMAL <<"first iteration automatically adjusted thresh="<<thresh
3308 0 : <<" ( "<<sigma<<"+adjustment factor:"<<adjFact<<")* rms )"<<LogIO::POST;
3309 : }
3310 : else {
3311 0 : thresh = 0.5*itsMax;
3312 : os << LogIO::NORMAL <<"first iteration automatically adjusted thresh="<<thresh
3313 0 : <<" (0.5*peak )"<<LogIO::POST;
3314 : }
3315 : }
3316 0 : if (adjFact != 0.0) {
3317 0 : os << LogIO::NORMAL <<"thresh="<<thresh<<" ( "<<sigma<<"+adjustment factor:"<<adjFact<<")* rms )"<<LogIO::POST;
3318 : }
3319 : else {
3320 0 : os << LogIO::NORMAL <<"thresh="<<thresh<<" ( "<<sigma<<"* rms )"<<LogIO::POST;
3321 : }
3322 : }
3323 :
3324 :
3325 0 : itsRms = maxRmsVal;
3326 :
3327 0 : if (thresh > maxMaxVal) {
3328 0 : os << LogIO::WARN <<" The threshold value,"<<thresh<<" for making a mask is greater than max value in the image. No new mask will be added by automask."<< LogIO::POST;
3329 0 : tempMask->set(0.0);
3330 : }
3331 : else {
3332 :
3333 : // apply threshold to rebinned image to generate a temp image mask
3334 : // first run pruning by limiting n masks (npix=1 as it is already binned)
3335 0 : std::shared_ptr<ImageInterface<Float> > dummyim = pruneRegions(tempRebinnedIm, thresh, nmask, 1);
3336 :
3337 0 : os << LogIO::DEBUG1<<" threshold applied ="<<thresh<<LogIO::POST;
3338 : //cerr<<"dummyim shape="<<dummyim.get()->shape()<<endl;
3339 : //cerr<<"temprebinned shape="<<tempRebinnedIm.shape()<<endl;
3340 : //
3341 : //LatticeExpr<Float> tempthresh( iif( abs(tempRebinnedIm) > thresh, 1.0, 0.0) );
3342 0 : LatticeExpr<Float> tempthresh( iif( abs( *(dummyim.get()) ) > thresh, 1.0, 0.0) );
3343 : //os << LogIO::DEBUG1<<" copying the threshold image....."<<LogIO::POST;
3344 0 : tempMask->copyData(tempthresh);
3345 : }
3346 0 : return std::shared_ptr<ImageInterface<Float> >(tempMask);
3347 : }
3348 :
3349 348 : std::shared_ptr<ImageInterface<Float> > SDMaskHandler::convolveMask(const ImageInterface<Float>& inmask, const GaussianBeam& beam)
3350 : {
3351 1044 : LogIO os( LogOrigin("SDMaskHandler","convolveMask",WHERE) );
3352 348 : TempImage<Float>* tempIm = new TempImage<Float>(inmask.shape(), inmask.coordinates(), memoryToUse() );
3353 348 : tempIm->copyData(inmask);
3354 696 : std::shared_ptr<casacore::ImageInterface<float> > tempIm2_ptr(tempIm);
3355 : //DEBUG will be removed
3356 :
3357 696 : Vector<Quantity> convbeam(3);
3358 348 : convbeam[0] = beam.getMajor();
3359 348 : convbeam[1] = beam.getMinor();
3360 348 : convbeam[2] = beam.getPA();
3361 348 : os << LogIO::DEBUG1<<"convolve with a gaussian: bmaj="<<convbeam[0]<<"bmin="<<convbeam[1]<<" pa="<<convbeam[2]<<LogIO::POST;
3362 696 : Record dammyRec=Record();
3363 : //String convimname("temp_convim");
3364 1044 : Image2DConvolver<Float> convolver(tempIm2_ptr, &dammyRec, String(""), String(""), True);
3365 348 : convolver.setKernel("GAUSSIAN", convbeam[0], convbeam[1], convbeam[2]);
3366 348 : convolver.setAxes(std::make_pair(0,1));
3367 348 : convolver.setScale(Double(-1.0));
3368 348 : convolver.setSuppressWarnings(True);
3369 348 : auto outmask = convolver.convolve();
3370 696 : return outmask;
3371 : }
3372 :
3373 0 : std::shared_ptr<ImageInterface<Float> > SDMaskHandler::convolveMask(const ImageInterface<Float>& inmask, Int nxpix, Int nypix)
3374 : {
3375 0 : LogIO os( LogOrigin("SDMaskHandler","convolveMask",WHERE) );
3376 0 : TempImage<Float>* tempIm = new TempImage<Float>(inmask.shape(), inmask.coordinates(), memoryToUse() );
3377 0 : tempIm->copyData(inmask);
3378 0 : std::shared_ptr<casacore::ImageInterface<float> > tempIm2_ptr(tempIm);
3379 : //DEBUG will be removed
3380 0 : os << LogIO::DEBUG1<<"convolve with "<<nxpix<<" pix x "<<nypix<<" pix gaussian"<<LogIO::POST;
3381 :
3382 0 : Vector<Quantity> convbeam(3);
3383 0 : convbeam[0] = Quantity(nxpix, "pix");
3384 0 : convbeam[1] = Quantity(nypix, "pix");
3385 0 : convbeam[2] = Quantity(0.0, "deg");
3386 0 : Record dammyRec=Record();
3387 : //String convimname("temp_convim");
3388 0 : Image2DConvolver<Float> convolver(tempIm2_ptr, &dammyRec, String(""), String(""), True);
3389 0 : convolver.setKernel("GAUSSIAN", convbeam[0], convbeam[1], convbeam[2]);
3390 0 : convolver.setAxes(std::make_pair(0,1));
3391 0 : convolver.setScale(Double(-1.0));
3392 0 : convolver.setSuppressWarnings(True);
3393 0 : auto outmask = convolver.convolve();
3394 0 : return outmask;
3395 : }
3396 :
3397 0 : std::shared_ptr<casacore::ImageInterface<Float> > SDMaskHandler::pruneRegions(const ImageInterface<Float>& image, Double& thresh, Int nmask, Int npix)
3398 : {
3399 0 : LogIO os( LogOrigin("SDMaskHandler", "pruneRegions",WHERE) );
3400 0 : Bool debug(False);
3401 :
3402 0 : IPosition fullimShape=image.shape();
3403 0 : TempImage<Float>* fullIm = new TempImage<Float>(TiledShape(fullimShape, image.niceCursorShape()), image.coordinates(), memoryToUse());
3404 :
3405 0 : if (nmask==0 && npix==0 ) {
3406 : //No-op
3407 0 : os<<LogIO::DEBUG1<<"Skip pruning of mask regions"<<LogIO::POST;
3408 0 : fullIm->copyData(image);
3409 0 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3410 : }
3411 0 : os <<LogIO::DEBUG1<< "pruneRegions with nmask="<<nmask<<", size="<<npix<<" is applied"<<LogIO::POST;
3412 :
3413 0 : IPosition shp = image.shape();
3414 : //cerr<<"shp = "<<shp<<endl;
3415 0 : IPosition blc(4,0);
3416 0 : IPosition trc = shp-1;
3417 0 : Slicer sl(blc,trc,Slicer::endIsLast);
3418 0 : AxesSpecifier aspec(False);
3419 : // decomposer can only work for 2 and 3-dim images so need
3420 : // some checks the case for stokes axis > 1
3421 : // following works if stokes axis dim = 1
3422 0 : SubImage<Float>* subIm = new SubImage<Float>(image, sl, aspec, True);
3423 0 : RegionManager regMan;
3424 0 : CoordinateSystem cSys=subIm->coordinates();
3425 0 : regMan.setcoordsys(cSys);
3426 : //String strReg = "box[["+String::toString(blc[0])+"pix,"+String::toString(blc[1])+"pix], ["+String::toString(shp[0])+"pix,"+String::toString(shp[1])+"pix]]";
3427 : //cerr<<"strReg="<<strReg<<endl;
3428 : //RegionTextList CRTFList(cSys, strReg, shp);
3429 : //Record regRec = CRTFList.regionAsRecord();
3430 : //std::shared_ptr<casacore::SubImage<Float> > subIm = SubImageFactory<Float>::createSubImageRW(image, regRec, "", &os, aspec, False, False);
3431 :
3432 : //cerr <<" subIm.shape="<<subIm->shape()<<endl;
3433 0 : IPosition subimShape=subIm->shape();
3434 0 : uInt ndim = subimShape.nelements();
3435 : //std::shared_ptr<casacore::ImageInterface<float> > tempIm_ptr(subIm);
3436 0 : TempImage<Float>* tempIm = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates(), memoryToUse() );
3437 : // to search for both positive and negative components
3438 0 : tempIm->copyData(LatticeExpr<Float> (abs(*subIm)));
3439 0 : os << LogIO::NORMAL2 <<"Finding regions by ImageDecomposer..."<<LogIO::POST;
3440 : //use ImageDecomposer
3441 0 : Matrix<Quantity> blctrcs;
3442 0 : ImageDecomposer<Float> id(*tempIm);
3443 0 : id.setDeblend(True);
3444 0 : os << LogIO::DEBUG1<< "Deblend threshold="<<thresh<<LogIO::POST;
3445 0 : id.setDeblendOptions(thresh, 3, 1, 2); //nContour=3
3446 : //id.setDeblendOptions(thresh, 3, 1, 1); //nContour=3, naxis=1
3447 0 : id.setFit(False);
3448 0 : os << LogIO::DEBUG1<<"Now calling decomposeImage()..."<<LogIO::POST;
3449 0 : id.decomposeImage();
3450 0 : if (debug)
3451 0 : id.printComponents();
3452 0 : uInt nRegion = id.numRegions();
3453 0 : os << "Found " << nRegion <<" regions"<<LogIO::POST;
3454 0 : Block<IPosition> blcs(nRegion);
3455 0 : Block<IPosition> trcs(nRegion);
3456 0 : id.boundRegions(blcs, trcs);
3457 : //os << "Get comp list.."<<LogIO::POST;
3458 0 : Matrix<Float> clmat=id.componentList();
3459 : //os << "Get peaks.."<<LogIO::POST;
3460 0 : Vector<Float> peaks = clmat.column(0);
3461 : //cerr<<"peaks="<<peaks<<endl;
3462 0 : os << LogIO::DEBUG1<< "Sorting by peak fluxes..."<<LogIO::POST;
3463 : // sort
3464 0 : Vector<uInt> sortedindx;
3465 0 : Sort sorter;
3466 : //cerr<<"Setup sortKey..."<<endl;
3467 0 : sorter.sortKey(peaks.data(),TpFloat,0, Sort::Descending);
3468 : //cerr<<"do sort..."<<endl;
3469 0 : sorter.sort(sortedindx, peaks.nelements());
3470 : //os << "Sorting py peak flux DONE..."<<LogIO::POST;
3471 0 : os<< LogIO::DEBUG1<<"sortedindx="<<sortedindx<<LogIO::POST;
3472 : // FOR DEBUGGING
3473 0 : for (uInt j = 0; j < blcs.nelements(); j++) {
3474 0 : os<<" j="<<j<<" blcs="<<blcs[j]<<" trcs="<<trcs[j]<<LogIO::POST;
3475 : }
3476 0 : Vector<Int> absRel(ndim, RegionType::Abs);
3477 0 : PtrBlock<const WCRegion *> wbox;
3478 0 : uInt iSubComp=0;
3479 0 : uInt removeByNMask=0;
3480 0 : uInt removeBySize=0;
3481 0 : for (uInt icomp=0; icomp < sortedindx.nelements(); icomp++) {
3482 0 : Bool removeit(False);
3483 0 : Vector<Quantum<Double> > qblc(ndim);
3484 0 : Vector<Quantum<Double> > qtrc(ndim);
3485 0 : Vector<Double> wblc(ndim);
3486 0 : Vector<Double> wtrc(ndim);
3487 0 : Vector<Double> pblc(ndim);
3488 0 : Vector<Double> ptrc(ndim);
3489 : // pixel blcs and trcs
3490 0 : for (uInt i=0; i < ndim; i++) {
3491 0 : pblc[i] = (Double) blcs[sortedindx[icomp]][i];
3492 0 : ptrc[i] = (Double) trcs[sortedindx[icomp]][i];
3493 : }
3494 : // get blcs and trcs in world coord.
3495 0 : cSys.toWorld(wblc,pblc);
3496 0 : cSys.toWorld(wtrc,ptrc);
3497 0 : for (uInt i=0; i < ndim; i++) {
3498 0 : qblc[i] = Quantum<Double>(wblc[i], cSys.worldAxisUnits()(cSys.pixelAxisToWorldAxis(i)));
3499 0 : qtrc[i] = Quantum<Double>(wtrc[i], cSys.worldAxisUnits()(cSys.pixelAxisToWorldAxis(i)));
3500 : }
3501 :
3502 0 : if (npix > 0) {
3503 : //npix = area
3504 : //os<<"pruning regions by size < "<<npix<<LogIO::POST;
3505 0 : Int xboxdim = ptrc[0] - pblc[0];
3506 0 : Int yboxdim = ptrc[1] - pblc[1];
3507 : //if (( xboxdim < npix || yboxdim < npix ) && xboxdim*yboxdim < npix*npix ) {
3508 0 : if ( xboxdim*yboxdim < npix ) {
3509 0 : removeit = True;
3510 0 : removeBySize++;
3511 : }
3512 : }
3513 0 : if (nmask > 0 && icomp >= (uInt)nmask ) {
3514 0 : removeit=True;
3515 0 : removeByNMask++;
3516 : }
3517 :
3518 0 : if (removeit) {
3519 0 : wbox.resize(iSubComp+1);
3520 0 : wbox[iSubComp]= new WCBox (qblc, qtrc, cSys, absRel);
3521 0 : iSubComp++;
3522 0 : os << LogIO::DEBUG1<<"*** Removed region: "<<icomp<<" pblc="<<pblc<<" ptrc="<<ptrc<<LogIO::POST;
3523 : }
3524 : else {
3525 0 : os << LogIO::DEBUG1<<"Saved region: "<<icomp<<" pblc="<<pblc<<" ptrc="<<ptrc<<LogIO::POST;
3526 : }
3527 : } // for icomp
3528 :
3529 : //cerr<<"iSubComp="<<iSubComp<<endl;
3530 : //cerr<<"wbox.nelements="<<wbox.nelements()<<endl;
3531 0 : if (iSubComp>0) {
3532 0 : ImageRegion* boxImageRegion=regMan.doUnion(wbox);
3533 : //cerr<<"regionToMask ..."<<endl;
3534 0 : tempIm->copyData(*subIm);
3535 0 : regionToMask(*tempIm,*boxImageRegion, Float(0.0));
3536 : //cerr<<"Done regionToMask..."<<endl;
3537 0 : os <<"pruneRegions removed "<<iSubComp<<" regions from the mask image"<<LogIO::POST;
3538 0 : os <<" (reasons: "<< removeBySize<<" by size "<<", "<<removeByNMask<<" by nmask)" <<LogIO::POST;
3539 0 : for (uInt k=0; k < wbox.nelements(); k++) {
3540 0 : delete wbox[k];
3541 : }
3542 : }
3543 : else {
3544 0 : os <<"No regions are removed by pruning" << LogIO::POST;
3545 : }
3546 : // Debug
3547 0 : if (debug) {
3548 0 : PagedImage<Float> debugPrunedIm(tempIm->shape(),tempIm->coordinates(),"prunedSub.Im");
3549 0 : debugPrunedIm.copyData(*tempIm);
3550 : }
3551 : //
3552 : // Inserting pruned result image to the input image
3553 0 : Array<Float> subimData;
3554 : //IPosition fullimShape=image.shape();
3555 : //TempImage<Float>* fullIm = new TempImage<Float>(TiledShape(fullimShape, image.niceCursorShape()), image.coordinates());
3556 0 : fullIm->set(0);
3557 0 : IPosition start(fullimShape.nelements(),0);
3558 0 : IPosition stride(fullimShape.nelements(),1);
3559 0 : if (ndim ==3) {
3560 0 : IPosition substart(3,0);
3561 0 : IPosition subshape(3,subimShape(0),subimShape(1),1);
3562 0 : IPosition substride(3,1,1,1);
3563 0 : uInt nchan=subimShape(2);
3564 : //cerr<<"shape tempIm ="<<tempIm->shape()<<endl;
3565 : //cerr<<"shape fullIm ="<<fullIm->shape()<<endl;
3566 0 : for (uInt ich=0; ich < nchan; ich++) {
3567 0 : substart(2) = ich;
3568 : //tempIm->getSlice(subimData,Slicer(substart,subend),True);
3569 0 : tempIm->getSlice(subimData,substart,subshape,substride,True);
3570 0 : start(3) = ich;
3571 0 : fullIm->putSlice(subimData,start,stride);
3572 : }
3573 : }
3574 0 : else if (ndim==2) {
3575 0 : subimData = tempIm->get();
3576 : //cerr<<"subimData shape="<<subimData.shape()<<endl;
3577 : //cerr<<"shape tempIm ="<<tempIm->shape()<<endl;
3578 0 : fullIm->putSlice(subimData,start,stride);
3579 : //cerr<<"shape fullIm ="<<fullIm->shape()<<endl;
3580 : }
3581 0 : delete subIm; subIm=0;
3582 0 : delete tempIm; tempIm=0;
3583 0 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3584 : }
3585 :
3586 :
3587 0 : std::shared_ptr<casacore::ImageInterface<Float> > SDMaskHandler::pruneRegions2(const ImageInterface<Float>& image, Double& thresh, Int nmask, Double prunesize)
3588 : {
3589 0 : LogIO os( LogOrigin("SDMaskHandler", "pruneRegions2",WHERE) );
3590 0 : Bool debug(False);
3591 :
3592 0 : IPosition fullimShape=image.shape();
3593 0 : TempImage<Float>* fullIm = new TempImage<Float>(TiledShape(fullimShape, image.niceCursorShape()), image.coordinates(), memoryToUse());
3594 0 : fullIm->set(0);
3595 :
3596 0 : if (nmask==0 && prunesize==0.0 ) {
3597 : //No-op
3598 0 : os<<LogIO::DEBUG1<<"Skip pruning of mask regions"<<LogIO::POST;
3599 0 : fullIm->copyData(image);
3600 0 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3601 : }
3602 0 : os <<LogIO::NORMAL<< "pruneRegions with nmask="<<nmask<<", size="<<prunesize<<" is applied"<<LogIO::POST;
3603 :
3604 0 : IPosition shp = image.shape();
3605 0 : IPosition trc = shp-1;
3606 0 : Int specaxis = CoordinateUtil::findSpectralAxis(image.coordinates());
3607 0 : uInt nchan = shp(specaxis);
3608 : // do a single channel plane at time
3609 0 : for (uInt ich = 0; ich < nchan; ich++) {
3610 0 : IPosition start(4, 0, 0, 0,ich);
3611 0 : IPosition length(4, shp(0),shp(1),shp(2),1);
3612 0 : Slicer sl(start, length);
3613 : //cerr<<"slicer sl ="<<sl<<endl;
3614 0 : AxesSpecifier aspec(False);
3615 : // decomposer can only work for 2 and 3-dim images so need
3616 : // some checks the case for stokes axis > 1
3617 : // following works if stokes axis dim = 1
3618 0 : SubImage<Float>* subIm = new SubImage<Float>(image, sl, aspec, True);
3619 0 : RegionManager regMan;
3620 0 : CoordinateSystem cSys=subIm->coordinates();
3621 0 : regMan.setcoordsys(cSys);
3622 0 : IPosition subimShape=subIm->shape();
3623 0 : uInt ndim = subimShape.nelements();
3624 0 : TempImage<Float>* tempIm = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates(), memoryToUse() );
3625 : // to search for both positive and negative components
3626 0 : tempIm->copyData(LatticeExpr<Float> (abs(*subIm)));
3627 : //use ImageDecomposer
3628 0 : Matrix<Quantity> blctrcs;
3629 0 : ImageDecomposer<Float> id(*tempIm);
3630 0 : id.setDeblend(True);
3631 0 : os << LogIO::DEBUG1<< "Deblend threshold="<<thresh<<LogIO::POST;
3632 0 : id.setDeblendOptions(thresh, 3, 1, 2); //nContour=3
3633 0 : id.setFit(False);
3634 0 : os << LogIO::DEBUG1<<"Now calling decomposeImage()..."<<LogIO::POST;
3635 0 : id.decomposeImage();
3636 0 : if (debug)
3637 0 : id.printComponents();
3638 0 : uInt nRegion = id.numRegions();
3639 0 : os << "Found " << nRegion <<" regions for channel plane "<<ich<<LogIO::POST;
3640 0 : if (nRegion) {
3641 0 : Block<IPosition> blcs(nRegion);
3642 0 : Block<IPosition> trcs(nRegion);
3643 0 : id.boundRegions(blcs, trcs);
3644 : //os << "Get comp list.."<<LogIO::POST;
3645 0 : Matrix<Float> clmat=id.componentList();
3646 : //os << "Get peaks.."<<LogIO::POST;
3647 0 : Vector<Float> peaks = clmat.column(0);
3648 : //cerr<<"peaks="<<peaks<<endl;
3649 0 : os << LogIO::DEBUG1<< "Sorting by peak fluxes..."<<LogIO::POST;
3650 : // sort
3651 0 : Vector<uInt> sortedindx;
3652 0 : Sort sorter;
3653 0 : sorter.sortKey(peaks.data(),TpFloat,0, Sort::Descending);
3654 0 : sorter.sort(sortedindx, peaks.nelements());
3655 : //os << "Sorting py peak flux DONE..."<<LogIO::POST;
3656 0 : os<< LogIO::DEBUG1<<"sortedindx="<<sortedindx<<LogIO::POST;
3657 : // FOR DEBUGGING
3658 0 : if (debug) {
3659 0 : for (uInt j = 0; j < blcs.nelements(); j++) {
3660 0 : os<<" j="<<j<<" blcs="<<blcs[j]<<" trcs="<<trcs[j]<<LogIO::POST;
3661 : }
3662 : }
3663 0 : Vector<Int> absRel(ndim, RegionType::Abs);
3664 0 : PtrBlock<const WCRegion *> wbox;
3665 0 : uInt iSubComp=0;
3666 0 : uInt removeByNMask=0;
3667 0 : uInt removeBySize=0;
3668 0 : for (uInt icomp=0; icomp < sortedindx.nelements(); icomp++) {
3669 0 : Bool removeit(False);
3670 0 : if (debug) {
3671 0 : cerr<<"sortedindx="<<sortedindx[icomp]<<" comp="<<clmat.row(sortedindx[icomp])<<endl;
3672 : }
3673 0 : Vector<Quantum<Double> > qblc(ndim);
3674 0 : Vector<Quantum<Double> > qtrc(ndim);
3675 0 : Vector<Double> wblc(ndim);
3676 0 : Vector<Double> wtrc(ndim);
3677 0 : Vector<Double> pblc(ndim);
3678 0 : Vector<Double> ptrc(ndim);
3679 : // pixel blcs and trcs
3680 0 : for (uInt i=0; i < ndim; i++) {
3681 0 : pblc[i] = (Double) blcs[sortedindx[icomp]][i];
3682 0 : ptrc[i] = (Double) trcs[sortedindx[icomp]][i];
3683 : }
3684 : // get blcs and trcs in world coord.
3685 0 : cSys.toWorld(wblc,pblc);
3686 0 : cSys.toWorld(wtrc,ptrc);
3687 0 : if (debug) {
3688 0 : cerr<<"cSys.workdAxisUnits=="<<cSys.worldAxisUnits()<<endl;
3689 0 : cerr<<"cSys increment = "<<cSys.increment()<<endl;
3690 : }
3691 0 : for (uInt i=0; i < ndim; i++) {
3692 0 : qblc[i] = Quantum<Double>(wblc[i], cSys.worldAxisUnits()(cSys.pixelAxisToWorldAxis(i)));
3693 0 : qtrc[i] = Quantum<Double>(wtrc[i], cSys.worldAxisUnits()(cSys.pixelAxisToWorldAxis(i)));
3694 : }
3695 0 : if (prunesize > 0.0) {
3696 : //npix = area
3697 0 : Int xboxdim = ptrc[0] - pblc[0] +1;
3698 0 : Int yboxdim = ptrc[1] - pblc[1] +1;
3699 : // get size from component : these seem to be in deg
3700 0 : Double ax1 = clmat.row(sortedindx[icomp])[3];
3701 0 : Double ax2 = ax1*clmat.row(sortedindx[icomp])[4];
3702 0 : Quantity qax1(ax1,"deg");
3703 0 : Quantity qax2(ax2,"deg");
3704 0 : Vector<Double> incVal = cSys.increment();
3705 0 : Vector<String> incUnit = cSys.worldAxisUnits();
3706 0 : Quantity qinc1(incVal[0],incUnit[0]);
3707 0 : Quantity qinc2(incVal[1],incUnit[1]);
3708 0 : Double pixAreaInRad = abs(qinc1.get(Unit("rad")).getValue() * qinc2.get(Unit("rad")).getValue());
3709 0 : Double regionInSR = C::pi * qax1.get(Unit("rad")).getValue() * qax2.get(Unit("rad")).getValue() / (4. * C::ln2);
3710 0 : Double regpix = regionInSR/pixAreaInRad;
3711 : //Double axpix1 = ceil(abs(qax1/(qinc1.convert(qax1),qinc1)).get().getValue());
3712 : //Double axpix2 = ceil(abs(qax2/(qinc2.convert(qax2),qinc2)).get().getValue());
3713 : //Int regpix = Int(C::pi * axpix1 * axpix2 /(4. * C::ln2));
3714 0 : if (debug) {
3715 0 : cerr<<"regpix="<<regpix<<" prunesize="<<prunesize<<" xboxdim="<<xboxdim<<" yboxdim="<<yboxdim<<endl;
3716 : }
3717 : // regpix seems to be a bit too small ..., try pruning by blc, trc ...
3718 0 : if ( xboxdim*yboxdim < Int(ceil(prunesize)) ) {
3719 : //if ( regpix < prunesize ) {
3720 0 : removeit = True;
3721 0 : removeBySize++;
3722 : }
3723 : }
3724 0 : if (nmask > 0 && icomp >= (uInt)nmask ) {
3725 0 : removeit=True;
3726 0 : removeByNMask++;
3727 : }
3728 :
3729 0 : if (removeit) {
3730 0 : wbox.resize(iSubComp+1);
3731 0 : wbox[iSubComp]= new WCBox (qblc, qtrc, cSys, absRel);
3732 0 : iSubComp++;
3733 0 : os << LogIO::DEBUG1<<"*** Removed region: "<<icomp<<" pblc="<<pblc<<" ptrc="<<ptrc<<LogIO::POST;
3734 : }
3735 : else {
3736 0 : os << LogIO::DEBUG1<<"Saved region: "<<icomp<<" pblc="<<pblc<<" ptrc="<<ptrc<<LogIO::POST;
3737 : }
3738 : } // for icomp
3739 :
3740 : //cerr<<"iSubComp="<<iSubComp<<endl;
3741 : //cerr<<"wbox.nelements="<<wbox.nelements()<<endl;
3742 0 : if (iSubComp>0) {
3743 0 : ImageRegion* boxImageRegion=regMan.doUnion(wbox);
3744 : //cerr<<"regionToMask ..."<<endl;
3745 0 : tempIm->copyData(*subIm);
3746 0 : regionToMask(*tempIm,*boxImageRegion, Float(0.0));
3747 : //cerr<<"Done regionToMask..."<<endl;
3748 0 : os <<"pruneRegions removed "<<iSubComp<<" regions from the mask image"<<LogIO::POST;
3749 0 : os <<" (reasons: "<< removeBySize<<" by size "<<", "<<removeByNMask<<" by nmask)" <<LogIO::POST;
3750 0 : for (uInt k=0; k < wbox.nelements(); k++) {
3751 0 : delete wbox[k];
3752 : }
3753 : }
3754 : else {
3755 0 : os <<"No regions are removed by pruning" << LogIO::POST;
3756 : }
3757 : // Debug
3758 0 : if (debug) {
3759 0 : PagedImage<Float> debugPrunedIm(tempIm->shape(),tempIm->coordinates(),"tmp-prunedSub.im");
3760 0 : debugPrunedIm.copyData(*tempIm);
3761 : }
3762 : //
3763 : // Inserting pruned result image to the input image
3764 0 : Array<Float> subimData;
3765 : //IPosition fullimShape=image.shape();
3766 : //TempImage<Float>* fullIm = new TempImage<Float>(TiledShape(fullimShape, image.niceCursorShape()), image.coordinates());
3767 :
3768 : //IPosition start(fullimShape.nelements(),0);
3769 : //IPosition stride(fullimShape.nelements(),1);
3770 : //IPosition substart(3,0);
3771 : //IPosition subshape(3,subimShape(0),subimShape(1),1);
3772 : //IPosition substride(3,1,1,1);
3773 : //tempIm->getSlice(subimData,Slicer(substart,subend),True);
3774 0 : tempIm->getSlice(subimData,IPosition(2,0), tempIm->shape(), IPosition(2,1,1));
3775 0 : fullIm->putSlice(subimData,start,IPosition(4,1,1,1,1));
3776 0 : delete tempIm; tempIm=0;
3777 0 : delete subIm; subIm=0;
3778 : }// if(nRegion) end
3779 : }
3780 0 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3781 : }
3782 :
3783 : //yet another pruneRegions - using connect component labelling with depth first search alogirthm ..
3784 343 : std::shared_ptr<casacore::ImageInterface<Float> > SDMaskHandler::YAPruneRegions(const ImageInterface<Float>& image, Vector<Bool>& chanflag, Vector<Bool>& allpruned, Vector<uInt>& nreg, Vector<uInt>& npruned, Double prunesize, Bool showchanlabel)
3785 : {
3786 1029 : LogIO os( LogOrigin("SDMaskHandler", "YAPruneRegions",WHERE) );
3787 343 : Timer timer;
3788 343 : Bool debug(False);
3789 343 : Bool recordPruned(False);
3790 343 : if (allpruned.nelements()>0) {
3791 232 : recordPruned=True;
3792 232 : allpruned.set(False);
3793 : }
3794 :
3795 686 : IPosition fullimShape=image.shape();
3796 343 : TempImage<Float>* fullIm = new TempImage<Float>(TiledShape(fullimShape, image.niceCursorShape()), image.coordinates(), memoryToUse());
3797 343 : fullIm->set(0);
3798 :
3799 343 : if (prunesize==0.0 ) {
3800 : //No-op
3801 0 : os<<LogIO::DEBUG1<<"Skip pruning of mask regions"<<LogIO::POST;
3802 0 : fullIm->copyData(image);
3803 0 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3804 : }
3805 343 : os <<LogIO::DEBUG1<< "pruneRegions with size="<<prunesize<<" is applied"<<LogIO::POST;
3806 :
3807 686 : IPosition shp = image.shape();
3808 343 : Int specaxis = CoordinateUtil::findSpectralAxis(image.coordinates());
3809 343 : uInt nchan = shp(specaxis);
3810 343 : nreg.resize(nchan);
3811 343 : npruned.resize(nchan);
3812 : // do a single channel plane at time
3813 : // - assumes standard CASA image axis ordering (ra,dec,stokes,chan)
3814 686 : for (uInt ich = 0; ich < nchan; ich++) {
3815 343 : if (!chanflag(ich)) {
3816 686 : IPosition start(4, 0, 0, 0,ich);
3817 686 : IPosition length(4, shp(0),shp(1),shp(2),1);
3818 686 : Slicer sl(start, length);
3819 : //cerr<<"ich="<<ich<<" slicer sl ="<<sl<<endl;
3820 686 : AxesSpecifier aspec(False);
3821 : // following works if stokes axis dim = 1
3822 343 : SubImage<Float>* subIm = new SubImage<Float>(image, sl, aspec, True);
3823 :
3824 686 : IPosition subimShape=subIm->shape();
3825 343 : TempImage<Float>* tempIm = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates(), memoryToUse() );
3826 : // to search for both positive and negative components
3827 343 : tempIm->copyData(LatticeExpr<Float> (abs(*subIm)));
3828 :
3829 343 : TempImage<Float>* blobMap = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates(), memoryToUse() );
3830 343 : blobMap->set(0);
3831 :
3832 : // connected componet labelling
3833 343 : os<<LogIO::DEBUG1<<"Calling labelRegions..."<<LogIO::POST;
3834 686 : Array<Float> tempImarr;
3835 343 : tempIm->get(tempImarr);
3836 343 : Float sumMaskVal=sum(tempImarr);
3837 343 : uInt removeBySize=0;
3838 343 : uInt nBlob=0;
3839 343 : os<<LogIO::DEBUG1<<" total pix of 1s="<< sumMaskVal <<LogIO::POST;
3840 343 : if ( sumMaskVal !=0.0 ) {
3841 309 : timer.mark();
3842 309 : labelRegions(*tempIm, *blobMap);
3843 309 : os<< LogIO::DEBUG1 << "Processing time for labelRegions: real "<< timer.real()<< "s ; user "<< timer.user() <<"s"<< LogIO::POST;
3844 618 : Array<Float> tempblobarr;
3845 309 : blobMap->get(tempblobarr);
3846 309 : os<<LogIO::DEBUG1<<" total pix of 1s="<< sum(tempblobarr) <<LogIO::POST;
3847 309 : os<<LogIO::DEBUG1<<"Calling findBlobSize..."<<LogIO::POST;
3848 : // get blobsizes (the vector contains each labeled region size (label # = ith element+1)
3849 : //timer.mark();
3850 618 : Vector<Float> blobsizes = findBlobSize(*blobMap);
3851 : //cerr<<"blobsizes="<<blobsizes<<endl;
3852 : //use ImageDecomposer
3853 : // book keeping of no of removed components`
3854 : //cerr<<"blobsizes.nelements()="<<blobsizes.nelements()<<endl;
3855 : //removing operations
3856 309 : nBlob = blobsizes.nelements();
3857 309 : if (blobsizes.nelements()) {
3858 309 : if (prunesize > 0.0) {
3859 622 : for (uInt icomp = 0; icomp < blobsizes.nelements(); ++icomp) {
3860 313 : if ( blobsizes[icomp] < prunesize ) {
3861 64 : Float blobid = Float(icomp+1);
3862 64 : removeBySize++;
3863 64 : tempIm->copyData( (LatticeExpr<Float>)( iif(*blobMap == blobid, 0.0, *tempIm ) ) );
3864 : }
3865 : }//for-loop
3866 : }
3867 : }
3868 : } // if-sumMaskVal!=0
3869 : // log reporting ...
3870 686 : String chanlabel("");
3871 343 : if (showchanlabel) {
3872 0 : chanlabel = "[C"+String::toString(ich)+"]";
3873 : }
3874 343 : if (removeBySize>0) {
3875 64 : os <<LogIO::DEBUG1<<chanlabel<<" pruneRegions removed "<<removeBySize<<" regions (out of "<<nBlob<<" ) from the mask image. "<<LogIO::POST;
3876 64 : if (recordPruned) {
3877 64 : if (removeBySize==nBlob) allpruned(ich) = True;
3878 : }
3879 : }
3880 : else {
3881 279 : if (sumMaskVal!=0.0) {
3882 245 : os <<LogIO::NORMAL<<chanlabel<<" No regions are removed in pruning process." << LogIO::POST;
3883 : }
3884 : else {
3885 34 : os <<LogIO::NORMAL<<chanlabel<<" No regions are found in this plane."<< LogIO::POST;
3886 : }
3887 : }
3888 343 : nreg[ich] = nBlob;
3889 343 : npruned[ich] = removeBySize;
3890 : // Debug
3891 343 : if (debug) {
3892 0 : PagedImage<Float> tempBlobMap(blobMap->shape(), blobMap->coordinates(), "tmp-Blob.map");
3893 0 : tempBlobMap.copyData(*blobMap);
3894 : }
3895 343 : Array<Float> subimData;
3896 343 : tempIm->getSlice(subimData,IPosition(2,0), tempIm->shape(), IPosition(2,1,1));
3897 343 : fullIm->putSlice(subimData,start,IPosition(4,1,1,1,1));
3898 343 : delete tempIm; tempIm=0;
3899 343 : delete subIm; subIm=0;
3900 343 : delete blobMap; blobMap=0;
3901 : } // if-skipmask
3902 : else {
3903 0 : nreg[ich] = 0;
3904 0 : npruned[ich] = 0;
3905 0 : if (showchanlabel) {
3906 0 : os<<LogIO::DEBUG1<<"Skipping chan "<<ich<<" from pruning"<<LogIO::POST;
3907 : }
3908 : else {
3909 0 : os<<LogIO::DEBUG1<<"Skipping this plane from pruning"<<LogIO::POST;
3910 : }
3911 :
3912 : }
3913 : }
3914 343 : return std::shared_ptr<ImageInterface<Float> >(fullIm);
3915 : }
3916 :
3917 :
3918 52 : Float SDMaskHandler::pixelBeamArea(const GaussianBeam& beam, const CoordinateSystem& csys)
3919 : {
3920 :
3921 104 : Quantity bmaj = beam.getMajor();
3922 104 : Quantity bmin = beam.getMinor();
3923 104 : Vector<Double> incVal = csys.increment();
3924 104 : Vector<String> incUnit = csys.worldAxisUnits();
3925 104 : Quantity qinc1(incVal[0],incUnit[0]);
3926 52 : Quantity qinc2(incVal[1],incUnit[1]);
3927 : // should in rad but make sure...
3928 52 : Double pixArea = abs(qinc1.get(Unit("rad")).getValue() * qinc2.get(Unit("rad")).getValue());
3929 52 : Double solidAngle = C::pi * bmaj.get(Unit("rad")).getValue() * bmin.get(Unit("rad")).getValue()/(4.* C::ln2);
3930 104 : return (Float) solidAngle/pixArea;
3931 :
3932 : }
3933 :
3934 20 : void SDMaskHandler::makePBMask(std::shared_ptr<SIImageStore> imstore, Float pblimit, Bool combinemask)
3935 : {
3936 60 : LogIO os( LogOrigin("SDMaskHandler","makePBMask",WHERE) );
3937 :
3938 20 : if( imstore->hasPB() ) // Projection algorithms will have this.
3939 : {
3940 23 : LatticeExpr<Float> themask;
3941 19 : if (combinemask && imstore->hasMask()) {
3942 7 : os<<"MAKE combined PB mask"<<LogIO::POST;
3943 7 : themask = LatticeExpr<Float> ( iif( (*(imstore->pb())) > pblimit, *(imstore->mask()), 0.0 ) );
3944 : }
3945 : else {
3946 12 : os<<"MAKE PB mask"<<LogIO::POST;
3947 24 : themask = LatticeExpr<Float> ( iif( (*(imstore->pb())) > pblimit , 1.0, 0.0 ) );
3948 : }
3949 15 : imstore->mask()->copyData( themask );
3950 : }
3951 : else // Calculate it here...
3952 : {
3953 : // Get antenna diameter
3954 : // Get frequency
3955 : // Assuming a Gaussian, construct a circle region at pblimit.
3956 :
3957 : // But for now...
3958 1 : throw(AipsError("Need PB/Sensitivity/Weight image before a PB-based mask can be made for "+imstore->getName()));
3959 :
3960 : }
3961 : // Also add option to just use the vpmanager or whatever centralized PB repository there will be (sometime in the distant future...).
3962 :
3963 15 : }// end of makePBMask
3964 :
3965 : //apply per channel plane threshold
3966 696 : void SDMaskHandler::makeMaskByPerChanThreshold(const ImageInterface<Float>& image, Vector<Bool>& chanflag, ImageInterface<Float>& mask, Vector<Float>& thresholds, Vector<Float>& masksizes)
3967 : {
3968 1392 : IPosition imshape = image.shape();
3969 :
3970 1392 : CoordinateSystem imcsys = image.coordinates();
3971 1392 : Vector<Int> diraxes = CoordinateUtil::findDirectionAxes(imcsys);
3972 696 : Int specaxis = CoordinateUtil::findSpectralAxis(imcsys);
3973 696 : uInt nchan = imshape (specaxis);
3974 696 : masksizes.resize(nchan);
3975 696 : if (nchan != thresholds.nelements()) {
3976 0 : throw(AipsError("Mismatch in the number of threshold values and the number of chan planes."));
3977 : }
3978 1392 : for (uInt ich=0; ich < nchan; ich++) {
3979 696 : if (!chanflag(ich)) {
3980 1392 : IPosition start(4, 0, 0, 0,ich);
3981 1392 : IPosition length(4, imshape(diraxes(0)),imshape(diraxes(1)),imshape(2),1);
3982 1392 : Slicer sl(start, length);
3983 :
3984 : // make a subImage for a channel slice
3985 1392 : AxesSpecifier aspec(False);
3986 1392 : SubImage<Float> chanImage(image, sl, aspec, true);
3987 696 : TempImage<Float>* tempChanImage = new TempImage<Float> (chanImage.shape(), chanImage.coordinates(), memoryToUse() );
3988 1392 : Array<Float> chanImageArr;
3989 696 : LatticeExpr<Float> chanMask;
3990 696 : if (thresholds(ich) < 0) {
3991 : //LatticeExpr<Float> chanMask(iif(chanImage < thresholds(ich),1.0, 0.0));
3992 5 : chanMask = LatticeExpr<Float> (iif(chanImage < thresholds(ich),1.0, 0.0));
3993 : }
3994 : else {
3995 : //LatticeExpr<Float> chanMask(iif(chanImage > thresholds(ich),1.0, 0.0));
3996 691 : chanMask = LatticeExpr<Float> (iif(chanImage > thresholds(ich),1.0, 0.0));
3997 : }
3998 696 : tempChanImage->copyData(chanMask);
3999 : //tempChanImage->getSlice(chanImageArr, IPosition(4,0), chanImage.shape(),IPosition(4,1,1,1,1));
4000 696 : tempChanImage->getSlice(chanImageArr, IPosition(2,0), chanImage.shape(),IPosition(2,1,1));
4001 696 : mask.putSlice(chanImageArr,start,IPosition(4,1,1,1,1));
4002 696 : masksizes[ich]=sum(chanImageArr);
4003 696 : delete tempChanImage; tempChanImage=0;
4004 : }
4005 : else {
4006 0 : masksizes[ich]=0;
4007 : // cerr<<"makeMaskByPerChanThresh: skipping chan="<<ich<<endl;
4008 : }
4009 : } // loop over chans
4010 696 : }
4011 :
4012 119 : void SDMaskHandler::binaryDilationCore(Lattice<Float>& inlattice,
4013 : Array<Float>& structure,
4014 : Lattice<Bool>& mask,
4015 : Array<Bool>& chanmask,
4016 : Lattice<Float>& outlattice)
4017 : {
4018 357 : LogIO os( LogOrigin("SDMaskHandler", "binaryDilation", WHERE) );
4019 : //
4020 : //IPosition cursorShape=inlattice.niceCursorShape();
4021 238 : IPosition inshape = inlattice.shape();
4022 119 : Int nx = inshape(0);
4023 119 : Int ny = inshape(1);
4024 : // assume here 3x3 structure elements (connectivity of either 1 or 2)
4025 238 : IPosition seshape = structure.shape();
4026 119 : Int se_nx =seshape(0);
4027 119 : Int se_ny =seshape(1);
4028 :
4029 119 : if (mask.shape()!=inshape) {
4030 0 : throw(AipsError("Incompartible mask shape. Need to be the same as the input image."));
4031 : }
4032 : // assume the origin of structure element is the center se(1,1)
4033 238 : IPosition cursorShape(4, nx, ny, 1, 1);
4034 238 : IPosition axisPath(4, 0, 1, 3, 2);
4035 : //cerr<<"cursorShape="<<cursorShape<<endl;
4036 : //cerr<<"structure="<<structure<<endl;
4037 238 : LatticeStepper tls(inlattice.shape(), cursorShape, axisPath);
4038 238 : RO_LatticeIterator<Float> li(inlattice, tls);
4039 238 : RO_LatticeIterator<Bool> mi(mask, tls);
4040 238 : LatticeIterator<Float> oli(outlattice,tls);
4041 : Int ich;
4042 238 : IPosition ipch(chanmask.shape().size(),0);
4043 :
4044 : // for debug
4045 : //Array<Float> initarr=inlattice.get();
4046 : //cerr<<"initarr sum pix="<<sum(initarr)<<endl;
4047 :
4048 238 : for (li.reset(), mi.reset(), oli.reset(), ich=0; !li.atEnd(); li++, mi++, oli++, ich++) {
4049 : //Array<Float> planeImage(li.cursor());
4050 238 : Array<Float> inMask(li.cursor());
4051 : //cerr<<"sum of inMask="<<sum(inMask)<<endl;
4052 238 : Array<Float> planeImage(inMask.shape());
4053 119 : planeImage.set(0);
4054 119 : planeImage=inMask;
4055 : //cerr<<"sum of planeImage before grow ="<<sum(planeImage)<<endl;
4056 238 : Array<Bool> planeMask(mi.cursor());
4057 119 : ipch(0)=ich;
4058 : // if masks are true do binary dilation...
4059 119 : if (ntrue(planeMask)>0 && chanmask(ipch)) {
4060 : //cerr<<"planeImage.shape()="<<planeImage.shape()<<endl;
4061 9687 : for (Int i=0; i < nx; i++) {
4062 1869600 : for (Int j=0; j < ny; j++) {
4063 1860000 : if (planeImage(IPosition(4,i,j,0,0))==1.0 && planeMask(IPosition(4,i,j,0,0)) ) {
4064 : //cerr<<"if planeImage ==1 i="<<i<<" j="<<j<<endl;
4065 : // Set the value for se(1,1)
4066 5780 : planeImage(IPosition(4,i,j,0,0)) = 2.0;
4067 23120 : for (Int ise=0; ise < se_nx; ise++) {
4068 69360 : for (Int jse = 0; jse < se_ny; jse++) {
4069 52020 : Int relx_se = ise - 1;
4070 52020 : Int rely_se = jse - 1;
4071 52020 : if (structure(IPosition(2,ise,jse)) && !(ise==1.0 && jse==1.0)) {
4072 : //cerr<<"structure("<<ise<<","<<jse<<")="<<structure(IPosition(2,ise,jse))<<endl;
4073 23120 : if ((i + relx_se) >= 0 && (i + relx_se) < nx &&
4074 23120 : (j + rely_se) >= 0 && (j + rely_se) < ny) {
4075 23120 : if (planeImage(IPosition(4,i+relx_se,j+rely_se,0,0))==0 ) {
4076 : // set to 2.0 to make distinction with the orignal 1.0 pixels
4077 106 : planeImage(IPosition(4, i+relx_se, j+rely_se,0,0))=2.0;
4078 : //cerr<<" i+relx_se="<<i+relx_se<<" j+rely_se="<<j+rely_se<<endl;
4079 : }
4080 : }
4081 : } // if(se(ise,jse)
4082 : } // if planeImage(i,j)=1.0
4083 : } // S.E. col loop
4084 : } // S.E. row loop
4085 : } // image col loop
4086 : } //inage row loop
4087 :
4088 :
4089 9687 : for (Int ii=0; ii < nx; ii++) {
4090 1869600 : for (Int jj=0; jj < ny; jj++) {
4091 1860000 : if (planeImage(IPosition(4,ii,jj,0,0))==2)
4092 5886 : planeImage(IPosition(4,ii,jj,0,0))=1;
4093 : }
4094 : }
4095 : } // if ntrure() ...
4096 119 : oli.woCursor() = planeImage;
4097 : }
4098 : //For debug
4099 : //Array<Float> afterinarr=inlattice.get();
4100 : //cerr<<"afaterinarr sum pix ="<<sum(afterinarr)<<endl;
4101 : //Array<Float> outarr = outlattice.get();
4102 : //cerr<<"outlattice sum pix ="<<sum(outarr)<<endl;
4103 119 : }
4104 :
4105 111 : void SDMaskHandler::binaryDilation(ImageInterface<Float>& inImage,
4106 : Array<Float>& structure,
4107 : Int niteration,
4108 : Lattice<Bool>& mask,
4109 : Array<Bool>& chanmask,
4110 : ImageInterface<Float>& outImage)
4111 : {
4112 333 : LogIO os( LogOrigin("SDMaskHandler", "binaryDilation", WHERE) );
4113 222 : ArrayLattice<Float> templattice(inImage.shape());
4114 111 : templattice.copyData(inImage);
4115 333 : TempImage<Float> diffTempImage(outImage.shape(), outImage.coordinates(), memoryToUse());
4116 111 : diffTempImage.set(1);
4117 : // initial grow mask
4118 111 : binaryDilationCore(inImage,structure,mask,chanmask,outImage);
4119 333 : LatticeExpr<Float> diffIm0( abs(templattice - outImage ) );
4120 :
4121 : // if the initial grow does not change mask (i.e. diffIm0 = 0)
4122 : // then it won't enter the while loop below.
4123 111 : diffTempImage.copyData(diffIm0);
4124 111 : Int iter = 1;
4125 119 : while (iter < niteration && !isEmptyMask(diffTempImage)) {
4126 8 : templattice.copyData(outImage);
4127 8 : binaryDilationCore(templattice,structure,mask,chanmask,outImage);
4128 16 : LatticeExpr<Float> diffIm( abs(templattice - outImage ) );
4129 8 : diffTempImage.copyData(diffIm);
4130 : /***
4131 : if (isEmptyMask(diffTempImage)) {
4132 : cerr<<"current iter"<<iter<<" diffim is 0 "<<endl;
4133 : }
4134 : else {
4135 : cerr<<"current iter"<<iter<<endl;
4136 : }
4137 : ***/
4138 8 : iter++;
4139 : }
4140 111 : os<<"grow iter done="<<iter<<LogIO::POST;
4141 111 : }
4142 :
4143 :
4144 0 : void SDMaskHandler::autoMaskWithinPB(std::shared_ptr<SIImageStore> imstore,
4145 : ImageInterface<Float>& posmask,
4146 : const Int iterdone,
4147 : Vector<Bool>& chanflag,
4148 : Record& robuststatsrec,
4149 : const String& alg,
4150 : const String& threshold,
4151 : const Float& fracofpeak,
4152 : const String& resolution,
4153 : const Float& resbybeam,
4154 : const Int nmask,
4155 : const Bool autoadjust,
4156 : const Float& sidelobethreshold,
4157 : const Float& noisethreshold,
4158 : const Float& lownoisethreshold,
4159 : const Float& negativethreshold,
4160 : const Float& cutthreshold,
4161 : const Float& smoothfactor,
4162 : const Float& minbeamfrac,
4163 : const Int growiterations,
4164 : const Bool dogrowprune,
4165 : const Float& minpercentchange,
4166 : const Bool verbose,
4167 : const Bool fastnoise,
4168 : const Bool isthresholdreached,
4169 : Float pblimit)
4170 : {
4171 0 : LogIO os( LogOrigin("SDMaskHandler","autoMaskWithinPB",WHERE) );
4172 : //Bool debug(False);
4173 :
4174 0 : os <<LogIO::DEBUG1<<"Calling autoMaskWithinPB .."<<LogIO::POST;
4175 : // changed to do automask ater pb mask is generated so automask do stats within pb mask
4176 0 : autoMask( imstore, posmask, iterdone, chanflag, robuststatsrec, alg, threshold, fracofpeak, resolution, resbybeam, nmask, autoadjust,
4177 : sidelobethreshold, noisethreshold, lownoisethreshold, negativethreshold, cutthreshold, smoothfactor,
4178 : minbeamfrac, growiterations, dogrowprune, minpercentchange, verbose, fastnoise, isthresholdreached, pblimit);
4179 :
4180 0 : if( imstore->hasPB() )
4181 : {
4182 0 : LatticeExpr<Float> themask( iif( (*(imstore->pb())) > pblimit , (*(imstore->mask())), 0.0 ) );
4183 0 : imstore->mask()->copyData( themask );
4184 :
4185 : /****
4186 : // apply pb mask as pixel mask for now. This will be converted to 1/0 image later
4187 : LatticeExpr<Bool> mask( iif(*(imstore->pb()) > pblimit, True, False));
4188 : os <<"calling MakeMask, hasPixelMask? "<<imstore->mask()->hasPixelMask()<<LogIO::POST;
4189 : os <<"calling MakeMask, hasRegion mask0? "<<imstore->mask()->hasRegion("mask0")<<LogIO::POST;
4190 : os <<"defaultMask "<<imstore->mask()->getDefaultMask()<<LogIO::POST;
4191 : //ImageRegion outreg=imstore->mask()->makeMask("mask0", False, True);
4192 : ImageRegion outreg=imstore.get()->mask()->makeMask("mask0", True, True);
4193 : LCRegion& outmask=outreg.asMask();
4194 : outmask.copyData(mask);
4195 : os <<"Before defineRegion"<<LogIO::POST;
4196 : imstore.get()->mask()->defineRegion("mask0", outreg, RegionHandler::Masks, True);
4197 : os <<"setDefMask"<<LogIO::POST;
4198 : imstore.get()->mask()->setDefaultMask("mask0");
4199 : imstore.get()->releaseImage(imstore.get()->mask());
4200 : if (debug) {
4201 : cerr<<"Make a copy"<<endl;
4202 : PagedImage<Float> temppb(imstore->mask()->shape(), imstore->mask()->coordinates(),"tempPB.Im");
4203 : temppb.copyData(*(imstore->mask()));
4204 : temppb.defineRegion("mask0", outreg, RegionHandler::Masks, True);
4205 : temppb.setDefaultMask("mask0");
4206 : }
4207 : ***/
4208 :
4209 : }
4210 : //autoMask( imstore, alg, threshold, fracofpeak, resolution, resbybeam, nmask);
4211 :
4212 : // else... same options as makePBMask (put it into a helper function)
4213 0 : }// end of autoMaskWithinPB
4214 :
4215 : //region labelling code
4216 0 : void SDMaskHandler::depthFirstSearch(Int x,
4217 : Int y,
4218 : Int cur_label,
4219 : Array<Float>& inlatarr,
4220 : Array<Float>& lablatarr)
4221 : {
4222 0 : Vector<Int> dx(4);
4223 0 : Vector<Int> dy(4);
4224 : // 4-direction connectivity
4225 0 : dx(0) = 1; dx(1)=0;dx(2)=-1;dx(3)=0;
4226 0 : dy(0) = 0; dy(1)=1;dy(2)=0;dy(3)=-1;
4227 : //IPosition inshape = inlat.shape();
4228 0 : IPosition inshape = inlatarr.shape();
4229 0 : Int nrow = inshape(0);
4230 0 : Int ncol = inshape(1);
4231 : // out of bound condition
4232 0 : if(x < 0 || x == nrow) return;
4233 0 : if(y < 0 || y == ncol) return;
4234 : //2d lattice is assumed
4235 0 : IPosition loc(2,x,y);
4236 : // already labelled or not value 1 pixel
4237 0 : if(lablatarr(loc) || !inlatarr(loc)) return;
4238 :
4239 : //cerr<<"cur_label="<<cur_label<<" loc="<<loc<<endl;
4240 :
4241 0 : lablatarr(loc) = Float(cur_label);
4242 : //lablat.putAt(Float(cur_label), loc);
4243 : //
4244 : //
4245 : //recursively check the neighbor
4246 0 : for (uInt inc = 0; inc < 4; ++inc)
4247 0 : depthFirstSearch(x + dx[inc], y + dy[inc], cur_label, inlatarr, lablatarr);
4248 : }
4249 :
4250 313 : void SDMaskHandler::depthFirstSearch2(Int x,
4251 : Int y,
4252 : Int cur_label,
4253 : Array<Float>& inlatarr,
4254 : Array<Float>& lablatarr)
4255 :
4256 : {
4257 313 : std::stack<IPosition> mystack;
4258 313 : IPosition inshape = inlatarr.shape();
4259 313 : Int nrow = inshape(0);
4260 313 : Int ncol = inshape(1);
4261 : // out of bound condition
4262 313 : if(x < 0 || x == nrow) return;
4263 313 : if(y < 0 || y == ncol) return;
4264 : //2d lattice is assumed
4265 313 : IPosition loc(2,x,y);
4266 : // if already visited or not mask region, return
4267 313 : if(lablatarr(loc) || !inlatarr(loc)) return;
4268 :
4269 626 : IPosition curloc;
4270 313 : mystack.push(loc);
4271 27995 : while (!mystack.empty()) {
4272 27682 : curloc = mystack.top();
4273 27682 : mystack.pop( );
4274 : //cerr<<"curloc="<<curloc<<" cur_label="<<cur_label<<endl;
4275 27682 : lablatarr(curloc) = Float(cur_label);
4276 55364 : Vector<IPosition> loclist = defineNeighbors(curloc, nrow, ncol);
4277 : //cerr<<"defineNeighbors done nelements="<<loclist.nelements()<<endl;
4278 27682 : if (loclist.nelements()) {
4279 138410 : for (uInt i=0; i < loclist.nelements(); ++i)
4280 : {
4281 110728 : if (inlatarr(loclist(i)) == 1 && lablatarr(loclist(i)) == 0 )
4282 : {
4283 27369 : mystack.push(loclist(i));
4284 : }
4285 : }
4286 : }
4287 : }
4288 : }
4289 :
4290 27682 : Vector<IPosition> SDMaskHandler::defineNeighbors(IPosition& pos, Int nx, Int ny)
4291 : {
4292 27682 : Vector<IPosition> neighbors(0);
4293 27682 : Int nelement=0;
4294 : // 4-direction connectivity
4295 55364 : Vector<Int> dx(4);
4296 55364 : Vector<Int> dy(4);
4297 27682 : dx(0) = 1; dx(1)=0;dx(2)=-1;dx(3)=0;
4298 27682 : dy(0) = 0; dy(1)=1;dy(2)=0;dy(3)=-1;
4299 138410 : for (uInt inc = 0; inc < 4; ++inc) {
4300 221456 : IPosition newpos(2,0);
4301 110728 : newpos(0) = pos(0)+dx(inc);
4302 110728 : newpos(1) = pos(1)+dy(inc);
4303 110728 : if (newpos(0) >= 0 && newpos(0) < nx && newpos(1) >= 0 && newpos(1) < ny)
4304 : {
4305 110728 : nelement++;
4306 110728 : neighbors.resize(nelement,True);
4307 110728 : neighbors(nelement-1)=newpos;
4308 : }
4309 : }
4310 55364 : return neighbors;
4311 : }
4312 :
4313 309 : void SDMaskHandler::labelRegions(Lattice<Float>& inlat, Lattice<Float>& lablat)
4314 : {
4315 309 : Int blobId = 0;
4316 618 : IPosition inshape = inlat.shape();
4317 309 : Int nrow = inshape(0);
4318 309 : Int ncol = inshape(1);
4319 618 : Array<Float> inlatarr;
4320 618 : Array<Float> lablatarr;
4321 309 : inlat.get(inlatarr);
4322 309 : lablat.get(lablatarr);
4323 : //cerr<<"IN labelRegions:: inlat.shape="<<inlat.shape()<<" lablat.shape="<<lablat.shape()<<" nrow="<<nrow<<" ncol="<<ncol<<endl;
4324 :
4325 309 : if ( sum(inlatarr) !=0.0 ) {
4326 34809 : for (Int i = 0; i < nrow; ++i)
4327 : {
4328 7084500 : for (Int j = 0; j < ncol; ++j)
4329 : {
4330 : // evaluating elements with lattice seems to be very slow...
4331 : // changed to use Arrarys
4332 : //if (!lablat(IPosition(2,i,j)) && inlat(IPosition(2,i,j) ) )
4333 7050000 : if (!lablatarr(IPosition(2,i,j)) && inlatarr(IPosition(2,i,j) ) )
4334 : //depthFirstSearch(i, j, ++blobId, inlatarr, lablatarr);
4335 : // Use non-recursive version
4336 313 : depthFirstSearch2(i, j, ++blobId, inlatarr, lablatarr);
4337 : }
4338 : }
4339 309 : lablat.put(lablatarr);
4340 : }
4341 : //cerr<<"done blobId="<<blobId<<endl;
4342 309 : }
4343 :
4344 309 : Vector<Float> SDMaskHandler::findBlobSize(Lattice<Float>& lablat)
4345 : {
4346 : // iterate through lablat (2D)
4347 : // find max label in lablat
4348 : // create groupsize list vector gsize(max-1)
4349 : // get val at each pixel in lablat (ival=lablat.get(loc)) and add 1 to gsize(ival-1)
4350 : // print each labelled comp's size...
4351 :
4352 927 : LogIO os( LogOrigin("SDMaskHandler","findBlobSize",WHERE) );
4353 618 : IPosition inshape = lablat.shape();
4354 309 : Int nrow = inshape(0);
4355 309 : Int ncol = inshape(1);
4356 : // getting max value via LatticeExprNode seems to be slower
4357 : //LatticeExprNode leMax=max(lablat);
4358 : //Float maxlab = leMax.getFloat();
4359 618 : Array<Float> lablatarr;
4360 309 : lablat.get(lablatarr);
4361 309 : Float maxlab = max(lablatarr);
4362 : //os<<LogIO::DEBUG1<<"maxlab="<<maxlab<<LogIO::POST;
4363 :
4364 309 : if (maxlab < 1.0) {
4365 0 : return Vector<Float>();
4366 : }
4367 618 : Vector<Float> blobsizes(Int(maxlab),0);
4368 34809 : for (Int i = 0; i < nrow; ++i)
4369 : {
4370 7084500 : for (Int j =0; j < ncol; ++j)
4371 : {
4372 : //IPosition loc(4, i, j, 0, 0);
4373 14100000 : IPosition loc(2, i, j);
4374 : //os<<LogIO::DEBUG1<<"i="<<i<<" j="<<j<<" labelat(loc)="<<lablat(loc)<<LogIO::POST;
4375 : //if (lablat(loc)) blobsizes[Int(lablat(loc))-1]+=1;
4376 7050000 : if (lablatarr(loc)) blobsizes[Int(lablatarr(loc))-1]+=1;
4377 : }
4378 : }
4379 :
4380 : //for debug
4381 622 : for (Int k = 0;k < maxlab; ++k)
4382 : {
4383 313 : os<<LogIO::DEBUG1<<"blobsizes["<<k<<"]="<<blobsizes[k]<<LogIO::POST;
4384 : }
4385 :
4386 309 : return blobsizes;
4387 : }
4388 :
4389 4 : void SDMaskHandler::printAutomaskSummary (const Array<Double>& rmss,
4390 : const Array<Double>& maxs,
4391 : const Array<Double>& mins,
4392 : const Array<Double>& mdns,
4393 : //const Vector<Float>& maskthreshold,
4394 : const Matrix<Float>& maskthreshold,
4395 : //const Vector<String>& thresholdtype,
4396 : const Matrix<String>& thresholdtype,
4397 : //const Vector<Bool>& chanflag,
4398 : const Vector<Bool>& chanflag,
4399 : //const Vector<Bool>& /* zeroChanMask */,
4400 : const Matrix<Bool>& /* zeroChanMask */,
4401 : //const Vector<uInt>& nreg,
4402 : const Matrix<uInt>& nreg,
4403 : //const Vector<uInt>& npruned,
4404 : const Matrix<uInt>& npruned,
4405 : //const Vector<uInt>& ngrowreg,
4406 : const Matrix<uInt>& ngrowreg,
4407 : //const Vector<uInt>& ngrowpruned,
4408 : const Matrix<uInt>& ngrowpruned,
4409 : //const Vector<Float>& negmaskpixs,
4410 : const Matrix<Float>& negmaskpixs,
4411 : const Record& miscSummaryInfo)
4412 :
4413 : {
4414 12 : LogIO os( LogOrigin("SDMaskHandler","printAutomaskSummary",WHERE) );
4415 :
4416 : // miscSummaryInfo currently contains sidelobe level and pruneregionsize
4417 : // but these won't be printed out here now (these are printed out in the beginning).
4418 : // or alll the arguments maybe packed into record...
4419 : Float sidelobelevel;
4420 4 : miscSummaryInfo.get("sidelobelevel", sidelobelevel);
4421 : Float prunesize;
4422 4 : miscSummaryInfo.get("pruneregionsize", prunesize);
4423 :
4424 4 : os << LogIO::NORMAL <<"========== automask summary ==========" << LogIO::POST;
4425 : os << LogIO::NORMAL <<"chan masking? median RMS"<<" "
4426 : <<"peak thresh_type thresh_value "
4427 4 : <<"N_reg N_pruned N_grow N_grow_pruned N_neg_pix"<<LogIO::POST;
4428 :
4429 : //use maskthreshold shape to find npol and nchan. Maskthreshold(npol, nchan)
4430 4 : Int npol = maskthreshold.nrow();
4431 4 : Int nchan = maskthreshold.ncolumn();
4432 4 : os << LogIO::DEBUG1 << "npol="<<npol<< " nchan="<<nchan<<LogIO::POST;
4433 8 : IPosition statshp = rmss.shape();
4434 : // Note: stats record collapse the axis with 1
4435 4 : os<<LogIO::DEBUG1 <<"rmss shape="<< String::toString(statshp) <<LogIO::POST;
4436 : // For the historical reason it is called chanidx but it is an index for a single plane stats
4437 8 : IPosition chanidx = statshp;
4438 : //uInt ndim = rmss.ndim();
4439 :
4440 : //Int nchan = maskthreshold.nelements();
4441 12 : for (uInt ich = 0; ich < (uInt) nchan; ich++) {
4442 8 : if (chanidx.nelements()==1 ) {
4443 2 : chanidx(0) = ich;
4444 : }
4445 6 : else if(chanidx.nelements()==2) {
4446 : // to include stats in all stokes in a single line
4447 6 : chanidx(1) = ich;
4448 6 : chanidx(0) =0;
4449 : }
4450 16 : Vector<Double> peaks(npol);
4451 40 : for (uInt ipol = 0; ipol < (uInt) npol; ipol++) {
4452 : //Double peak = abs(maxs(chanidx))> abs( mins(chanidx))? maxs(chanidx): mins(chanidx);
4453 32 : if (npol!=1) {
4454 32 : chanidx(0) = ipol;
4455 : }
4456 32 : peaks(ipol) = abs(maxs(chanidx))> abs( mins(chanidx))? maxs(chanidx): mins(chanidx);
4457 : //os << LogIO::DEBUG1 << "chanidx="<<chanidx<< " peaks("<<ipol<<")="<<peaks(ipol)<<LogIO::POST;
4458 :
4459 : }
4460 16 : String peakStr = npol==1? String::toString(peaks(0)):String::toString(peaks);
4461 : // only tested for single pol (normally stokes I)
4462 16 : String domasking = chanflag(ich)==0? "T":"F";
4463 : //String domasking = zeroChanMask[ich]==1? "F":"T";
4464 16 : String mdnsStr, rmssStr, maskthresholdStr;
4465 16 : String Nreg, Npruned, Ngrowreg, NgrowPruned, Nnegpix;
4466 16 : String NAstr("--");
4467 :
4468 : //if masking is skipped median, rms, thresholdvalue are
4469 : //set to ---
4470 8 : if (domasking=="F") {
4471 0 : mdnsStr=NAstr;
4472 0 : rmssStr=NAstr;
4473 0 : maskthresholdStr=NAstr;
4474 : }
4475 : else {
4476 : //mdnsStr=String::toString(mdns(chanidx));
4477 : //rmssStr=String::toString(rmss(chanidx));
4478 : //maskthresholdStr=String::toString(maskthreshold[ich]);
4479 : //reset pol axis of chanidx
4480 16 : IPosition trc = chanidx;
4481 8 : if (npol > 1) {
4482 8 : chanidx(0)=0;
4483 8 : trc(0)=npol-1;
4484 : }
4485 :
4486 8 : os<<LogIO::DEBUG1<<" mdns.shape="<<mdns.shape()<<"chanidx="<<chanidx<<" trc="<<trc<<LogIO::POST;
4487 8 : if (chanidx.nelements()==1) {
4488 2 : if (nchan==1) {
4489 2 : mdnsStr=String::toString(mdns);
4490 2 : rmssStr=String::toString(rmss);
4491 : }
4492 : else {
4493 0 : mdnsStr=String::toString(mdns(chanidx));
4494 0 : rmssStr=String::toString(rmss(chanidx));
4495 : }
4496 : }
4497 : else {
4498 12 : Matrix<Double> mdnsChan=mdns(chanidx, trc);
4499 6 : Matrix<Double> rmssChan=rmss(chanidx, trc);
4500 : //mdnsStr=String::toString(mdns(chanidx,trc));
4501 : //rmssStr=String::toString(rmss(chanidx,trc));
4502 6 : mdnsStr=String::toString(mdnsChan.column(0));
4503 6 : rmssStr=String::toString(rmssChan.column(0));
4504 : }
4505 16 : Vector<Float> maskthreshvec = maskthreshold.column(ich);
4506 8 : if (maskthreshvec.nelements()==1) {
4507 0 : maskthresholdStr=String::toString(maskthreshvec(0));
4508 : }
4509 : else {
4510 8 : maskthresholdStr=String::toString(maskthreshvec);
4511 : }
4512 : }
4513 :
4514 8 : if (!nreg.nelements()) {
4515 0 : Nreg = NAstr;
4516 : }
4517 : else {
4518 : //Nreg = String::toString(nreg[ich]);
4519 16 : Vector<uInt> nregvec = nreg.column(ich);
4520 8 : if (nregvec.nelements()==1) {
4521 0 : Nreg = String::toString(nregvec(0));
4522 : }
4523 : else {
4524 8 : Nreg = String::toString(nregvec);
4525 : }
4526 : }
4527 8 : if (!npruned.nelements()) {
4528 0 : Npruned = NAstr;
4529 : }
4530 : else {
4531 : //Npruned = String::toString(npruned[ich]);
4532 16 : Vector<uInt> nprunedvec = npruned.column(ich);
4533 8 : if (nprunedvec.nelements()==1) {
4534 0 : Npruned = String::toString(nprunedvec(0));
4535 : }
4536 : else {
4537 8 : Npruned = String::toString(nprunedvec);
4538 : }
4539 : }
4540 8 : if (!ngrowreg.nelements()) {
4541 0 : Ngrowreg = NAstr;
4542 : }
4543 : else {
4544 : //Ngrowreg = String::toString(ngrowreg[ich]);
4545 16 : Vector<uInt> ngrowregvec = ngrowreg.column(ich);
4546 8 : if (ngrowregvec.nelements()==1) {
4547 0 : Ngrowreg = String::toString(ngrowregvec(0));
4548 : }
4549 : else {
4550 8 : Ngrowreg = String::toString(ngrowreg.column(ich));
4551 : }
4552 : }
4553 8 : if (!ngrowpruned.nelements()) {
4554 0 : NgrowPruned = NAstr;
4555 : }
4556 : else {
4557 : //NgrowPruned = String::toString(ngrowpruned[ich]);
4558 16 : Vector<uInt> ngrowprunedvec=ngrowpruned.column(ich);
4559 8 : if (ngrowprunedvec.nelements()==1) {
4560 0 : NgrowPruned = String::toString(ngrowprunedvec(0));
4561 : }
4562 : else {
4563 8 : NgrowPruned = String::toString(ngrowprunedvec);
4564 : }
4565 : }
4566 8 : if (!negmaskpixs.nelements()) {
4567 0 : Nnegpix = NAstr;
4568 : }
4569 : else {
4570 : //Nnegpix = String::toString(negmaskpixs[ich]);
4571 16 : Vector<Float> negmaskpixsvec=negmaskpixs.column(ich);
4572 8 : if (negmaskpixsvec.nelements()==1) {
4573 0 : Nnegpix = String::toString(negmaskpixsvec(0));
4574 : }
4575 : else {
4576 8 : Nnegpix = String::toString(negmaskpixsvec);
4577 : }
4578 : }
4579 16 : Vector<String> threshtypevec=thresholdtype.column(ich);
4580 16 : String threshtypeStr;
4581 8 : if (threshtypevec.nelements()==1) {
4582 0 : threshtypeStr=threshtypevec(0);
4583 : }
4584 : else {
4585 8 : threshtypeStr=String::toString(threshtypevec);
4586 : }
4587 :
4588 : os << LogIO::NORMAL << "[C" << ich << "] "
4589 : << domasking << " "
4590 : //<< mdns(chanidx) << " "
4591 : //<< rmss(chanidx) << " "
4592 : << mdnsStr << " "
4593 : << rmssStr << " "
4594 : << peakStr << " "
4595 : << threshtypeStr << " "
4596 : //<< maskthreshold[ich] << " "
4597 : << maskthresholdStr << " "
4598 : << Nreg << " "
4599 : << Npruned << " "
4600 : << Ngrowreg << " "
4601 : << NgrowPruned << " "
4602 : << Nnegpix
4603 8 : << LogIO::POST;
4604 : }
4605 4 : os << LogIO::NORMAL <<"========== END of automask summary ==========" << LogIO::POST;
4606 4 : }
4607 :
4608 143 : void SDMaskHandler::setPBMaskLevel (const Float pbmasklevel) {
4609 143 : itsPBMaskLevel = pbmasklevel;
4610 143 : }
4611 :
4612 0 : Float SDMaskHandler::getPBMaskLevel() {
4613 0 : return itsPBMaskLevel;
4614 : }
4615 :
4616 : } //# NAMESPACE CASA - END
|
