Line data Source code
1 : //# SDMaskHandler.h: Definition for SDMaskHandler
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
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 adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_SDMASKHANDLER_H
30 : #define SYNTHESIS_SDMASKHANDLER_H
31 :
32 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/images/Images/ImageInterface.h>
35 : #include <casacore/images/Images/PagedImage.h>
36 : #include <casacore/images/Images/TempImage.h>
37 : #include <casacore/casa/Logging/LogMessage.h>
38 : #include <casacore/casa/Logging/LogSink.h>
39 :
40 : #include<synthesis/ImagerObjects/SIImageStore.h>
41 : #include<synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
42 :
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 :
45 : class SDMaskHandler
46 : {
47 : public:
48 :
49 : // Empty constructor
50 : SDMaskHandler();
51 : ~SDMaskHandler();
52 :
53 : void resetMask(std::shared_ptr<SIImageStore> imstore);
54 :
55 : void fillMask(std::shared_ptr<SIImageStore> imstore, casacore::Vector<casacore::String> maskStrings);
56 : void fillMask(std::shared_ptr<SIImageStore> imstore, casacore::String maskString);
57 :
58 : // Collection of methods translate mask description (text, record, threshold, etc) to
59 : // mask image where the region(s) of interest are represented by the value (default = 1.0)
60 : // and the rest of the image is set to 0.0.
61 : //
62 : //void makeMask();
63 : // Create a mask image with maskName from tempim with a threshold applied to the pixel intensity
64 : std::shared_ptr<casacore::ImageInterface<casacore::Float> > makeMask(const casacore::String& maskName, const casacore::Quantity threshold, casacore::ImageInterface<casacore::Float>& tempim);
65 :
66 : // Make a mask image from casacore::Record, casacore::Matrix of (nboxes,4) where each row contains [blc_x,blc_y, trc_x,trc_y],
67 : // and casacore::Matrix of (ncircles, 3) with the specified 'value'. Each row of circles are radius, x_center, y_center in pixels.
68 : // blctrcs and circles are applied to all the spectral and stokes while regions specified in record can specify selections
69 : // in spectral and stokes axes.
70 : //static casacore::Bool regionToImageMask(const casacore::String& maskimage, casacore::Record* regionRec, casacore::Matrix<casacore::Quantity> & blctrcs,
71 : static casacore::Bool regionToImageMask(casacore::ImageInterface<casacore::Float>& maskImage, const casacore::Record* regionRec, const casacore::Matrix<casacore::Quantity> & blctrcs,
72 0 : const casacore::Matrix<casacore::Float>& circles, const casacore::Float& value=1.0);
73 :
74 : // Convert boxes defined with blcs and trcs to ImageRegion
75 : static void boxRegionToImageRegion(const casacore::ImageInterface<casacore::Float>& maskImage, const casacore::Matrix<casacore::Quantity>& blctrcs, casacore::ImageRegion*& boxImageRegions);
76 : // Convert circles (in pixels) to ImageRegion
77 : static void circleRegionToImageRegion(const casacore::ImageInterface<casacore::Float>& maskImage, const casacore::Matrix<casacore::Float>& circles, casacore::ImageRegion*& circleImageRegions);
78 : // Convert region defined by record to Imageregion
79 : static void recordRegionToImageRegion(const casacore::Record* imageRegRec, casacore::ImageRegion*& imageRegion );
80 : // Convert casacore::ImageRegion to a mask image with the value
81 : static casacore::Bool regionToMask(casacore::ImageInterface<casacore::Float>& maskImage, casacore::ImageRegion& imageregion, const casacore::Float& value);
82 : // Read CRTF format text or the text file contains CRTF definitions and convert it to a ImageRegion
83 : static void regionTextToImageRegion(const casacore::String& text, const casacore::ImageInterface<casacore::Float>& regionImage, casacore::ImageRegion*& imageRegion);
84 :
85 : // merge mask images to outimage
86 : void copyAllMasks(const casacore::Vector< std::shared_ptr<casacore::ImageInterface<casacore::Float> > > inImageMasks, casacore::ImageInterface<casacore::Float>& outImageMask);
87 : // copy and regrid a mask image to outimage
88 : void copyMask(const casacore::ImageInterface<casacore::Float>& inimage, casacore::ImageInterface<casacore::Float>& outimage);
89 : // expand smaller chan mask image to larger one. - currently only works for a single channel (continuum) input mask
90 : void expandMask(const casacore::ImageInterface<casacore::Float>& inImageMask, casacore::ImageInterface<casacore::Float>& outImageMask);
91 : // convert internal mask to imageRegion
92 : void InMaskToImageRegion(const casacore::ImageInterface<casacore::Float>& inimage );
93 :
94 : int makeInteractiveMask(std::shared_ptr<SIImageStore>& imstore,
95 : casacore::Int& niter, casacore::Int& cycleniter,
96 : casacore::String& threshold, casacore::String& cyclethreshold);
97 :
98 : // Return a reference to an imageinterface for the mask.
99 : void makeAutoMask(std::shared_ptr<SIImageStore> imstore);
100 : // Top level autoMask interface
101 : // Different automasking algorithm can be choosen by specifying a specific alogrithm name in alg.
102 : // Not all arguments are used for a choosen alogrithm.
103 : //
104 : //
105 : // @param[in,out] imstore SIImageStore
106 : // @param[in,out] positive only mask
107 : // @param[in] current iteration number completed
108 : // @param[in, out] channel flag
109 : // @param[in, out] robust image statistics
110 : // @param[in] alg autoboxing alogrithm name (currently recongnized names are 'one-box', 'thresh','thresh2', and 'multithresh')
111 : // @param[in] threshold threshold with a unit
112 : // @param[in] fracpeak Fraction of the peak
113 : // @param[in] resolution Resolution with a unit
114 : // @param[in] resbybeam Multiplier of a beam to specify the resolution
115 : // @param[in] nmask Number of masks to use
116 : // @param[in] autoadjust Use autoadjustment of mask threshold (for 'thresh')
117 : // @param[in] sidlobethreshold Threshold factor in a multiplier of the sidelobe level
118 : // @param[in] noisethreshold Threshold factor in a multiplier of the rms noise
119 : // @param[in] lownoisethreshold Threshold factor in a multiplier of the rms noise used in the binary dilation
120 : // @param[in] negativethreshold Threshold factor in a multiplier of the rms noise used to set threshold for negative features
121 : // @param[in] cutthreshold Cut threshold factor for adjust a mask after smoothing of the mask
122 : // @param[in] smoothfactor Smoothing factor (multiplier of the beam)
123 : // @param[in] minbeamfrac Percent change in mask size to trigger a new automask creation for 'noise'-based threshold
124 : // @param[in] growiterations (maximum) number of binary dilation iteartions to grow the mask
125 : // @param[in] dogrowprune Toggle to do or skip pruning on the grow mask
126 : // @param[in] minpercentchange Mininum percentage change in mask to stop updating mask
127 : // @param[in] verbose Controls automask related logging messages
128 : // @param[in] isthresholdreached Check if cyclethreshold reached threshold
129 : // @param[in] fastnoise Toggle to turn on and off fast (but less robust) noise calculation
130 : // @param[in] pblimit Primary beam cut off level
131 : //
132 : void autoMask(std::shared_ptr<SIImageStore> imstore,
133 : casacore::ImageInterface<casacore::Float>& posmask,
134 : const casacore::Int iterdone,
135 : casacore::Vector<casacore::Bool>& chanflag,
136 : casacore::Record& robuststatsrec,
137 : const casacore::String& alg="",
138 : const casacore::String& threshold="",
139 : const casacore::Float& fracpeak=0.0,
140 : const casacore::String& resolution="",
141 : const casacore::Float& resbybeam=0.0,
142 : const casacore::Int nmask=0,
143 : const casacore::Bool autoadjust=false,
144 : const casacore::Float& sidelobethreshold=0.0,
145 : const casacore::Float& noisethreshold=0.0,
146 : const casacore::Float& lownoisethreshold=0.0,
147 : const casacore::Float& negativethreshold=0.0,
148 : const casacore::Float& cutthreshold=0.0,
149 : const casacore::Float& smoothfactor=0.0,
150 : const casacore::Float& minbeamfrac=0.0,
151 : const casacore::Int growiterations=0,
152 : const casacore::Bool dogrowprune=true,
153 : const casacore::Float& minpercentchange=0.0,
154 : const casacore::Bool verbose=false,
155 : const casacore::Bool fastnoise=false,
156 : const casacore::Bool isthresholdreached=false,
157 : casacore::Float pblimit=0.0);
158 :
159 : // automask by threshold with binning before applying it
160 : void autoMaskByThreshold (casacore::ImageInterface<casacore::Float>& mask,
161 : const casacore::ImageInterface<casacore::Float>& res,
162 : const casacore::ImageInterface<casacore::Float>& psf,
163 : const casacore::Quantity& resolution,
164 : const casacore::Float& resbybeam,
165 : const casacore::Quantity& qthreshold,
166 : const casacore::Float& fracofpeak,
167 : const casacore::Record& theStatsffff,
168 : const casacore::Float& sigma=3.0,
169 : const casacore::Int nmask=0,
170 : const casacore::Bool autoadjust=casacore::False);
171 :
172 : // automask by threshold : no binning version
173 : void autoMaskByThreshold2 (casacore::ImageInterface<casacore::Float>& mask,
174 : const casacore::ImageInterface<casacore::Float>& res,
175 : const casacore::ImageInterface<casacore::Float>& psf,
176 : const casacore::Quantity& resolution,
177 : const casacore::Float& resbybeam,
178 : const casacore::Quantity& qthreshold,
179 : const casacore::Float& fracofpeak,
180 : const casacore::Record& theStats,
181 : const casacore::Float& sigma=3.0,
182 : const casacore::Int nmask=0);
183 :
184 : // implementation of Amanda's automasking algorithm using multiple thresholds
185 : void autoMaskByMultiThreshold(casacore::ImageInterface<float>& mask,
186 : casacore::ImageInterface<casacore::Float>& posmask,
187 : const casacore::ImageInterface<casacore::Float>& res,
188 : const casacore::ImageInterface<casacore::Float>& psf,
189 : const casacore::Record& stats,
190 : const casacore::Record& newstats,
191 : const casacore::Int iterdone,
192 : casacore::Vector<casacore::Bool>& chanFlag,
193 : const casacore::Float& maskPercentChange=0.0,
194 : const casacore::Float& sidelobeLevel=0.0,
195 : const casacore::Float& sidelobeThresholdFactor=3.0,
196 : const casacore::Float& noiseThresholdFactor=3.0,
197 : const casacore::Float& lowNoiseThresholdFactor=2.0,
198 : const casacore::Float& negativeThresholdFactor=0.0,
199 : const casacore::Float& cutThreshold=0.01,
200 : const casacore::Float& smoothFactor=1.0,
201 : const casacore::Float& minBeamFrac=-1.0,
202 : const casacore::Int growIterations=100,
203 : const casacore::Bool dogrowprune=true,
204 : const casacore::Bool verbose=false,
205 : const casacore::Bool isthresholdreached=false);
206 :
207 :
208 : std::shared_ptr<casacore::ImageInterface<float> > makeMaskFromBinnedImage (
209 : const casacore::ImageInterface<casacore::Float>& image,
210 : const casacore::Int nx,
211 : const casacore::Int ny,
212 : const casacore::Float& fracofpeak,
213 : const casacore::Float& sigma,
214 : const casacore::Int nmask,
215 : const casacore::Bool autoadjust,
216 : casacore::Double thresh=0.0);
217 :
218 : // Convolve mask image with nx pixel by ny pixel
219 : std::shared_ptr<casacore::ImageInterface<float> > convolveMask(const casacore::ImageInterface<casacore::Float>& inmask,
220 : casacore::Int nxpix, casacore::Int nypix);
221 :
222 : // Convolve mask image by a gaussian
223 : std::shared_ptr<casacore::ImageInterface<float> > convolveMask(const casacore::ImageInterface<casacore::Float>& inmask,
224 : const casacore::GaussianBeam& beam);
225 :
226 : // Prune the mask regions found
227 : std::shared_ptr<casacore::ImageInterface<float> > pruneRegions(const casacore::ImageInterface<casacore::Float>& image,
228 : casacore::Double& thresh,
229 : casacore::Int nmask=0,
230 : casacore::Int npix=0);
231 : // Prune the mask regions per spectral plane
232 : std::shared_ptr<casacore::ImageInterface<float> > pruneRegions2(const casacore::ImageInterface<casacore::Float>& image,
233 : casacore::Double& thresh,
234 : casacore::Int nmask=0,
235 : casacore::Double prunesize=0.0);
236 :
237 : // Yet another Prune the mask regions per spectral plane
238 : std::shared_ptr<casacore::ImageInterface<float> > YAPruneRegions(const casacore::ImageInterface<casacore::Float>& image,
239 : casacore::Vector<casacore::Bool>& chanflag,
240 : casacore::Vector<casacore::Bool>& allpruned,
241 : casacore::Vector<casacore::uInt>& nreg,
242 : casacore::Vector<casacore::uInt>& npruned,
243 : casacore::Double prunesize=0.0,
244 : casacore::Bool showchanlabel=true);
245 :
246 : // create a mask image (1/0 image) applying a different threshold for each channel plane
247 : void makeMaskByPerChanThreshold(const casacore::ImageInterface<casacore::Float>& image,
248 : casacore::Vector<casacore::Bool>& chanflag,
249 : casacore::ImageInterface<casacore::Float>& mask,
250 : casacore::Vector<casacore::Float>& thresholds,
251 : casacore::Vector<casacore::Float>& masksizes);
252 :
253 : // A core method for binary dilation of the input lattice
254 : void binaryDilationCore(casacore::Lattice<casacore::Float>& inlattice,
255 : casacore::Array<casacore::Float>& structure,
256 : casacore::Lattice<casacore::Bool>& mask,
257 : casacore::Array<casacore::Bool>& chanmask,
258 : casacore::Lattice<casacore::Float>& outlattice);
259 :
260 : // Multiple Binary dilation application of an image with a constraint mask and channel plane based flags
261 : void binaryDilation(casacore::ImageInterface<casacore::Float>& inImage,
262 : casacore::Array<casacore::Float>& structure,
263 : casacore::Int niteration,
264 : casacore::Lattice<casacore::Bool>& mask,
265 : casacore::Array<casacore::Bool>& chanmask,
266 : casacore::ImageInterface<casacore::Float>& outImage);
267 :
268 : // return beam area in pixel unit
269 : casacore::Float pixelBeamArea(const casacore::GaussianBeam& beam, const casacore::CoordinateSystem& csys);
270 :
271 : // Create a mask image applying PB level
272 : // @param[in, out] imstore SIImageStore
273 : // @param[in] pblimit Primary beam cut off level
274 : // @param[in] pblimit Primary beam cut off level
275 : void makePBMask(std::shared_ptr<SIImageStore> imstore, casacore::Float pblimit=0.1, casacore::Bool combinemask=false);
276 :
277 : void autoMaskWithinPB(std::shared_ptr<SIImageStore> imstore,
278 : casacore::ImageInterface<casacore::Float>& posmask,
279 : const casacore::Int iterdone,
280 : casacore::Vector<casacore::Bool>& chanflag,
281 : casacore::Record& robuststatsrec,
282 : const casacore::String& alg="",
283 : const casacore::String& threshold="",
284 : const casacore::Float& fracpeak=0.0,
285 : const casacore::String& resolution="",
286 : const casacore::Float& resbybeam=0.0,
287 : const casacore::Int nmask=0,
288 : const casacore::Bool autoadjust=false,
289 : const casacore::Float& sidelobethreshold=0.0,
290 : const casacore::Float& noisethreshold=0.0,
291 : const casacore::Float& lownoisethreshold=0.0,
292 : const casacore::Float& negativethreshold=0.0,
293 : const casacore::Float& cutthreshold=0.0,
294 : const casacore::Float& smoothfactor=0.0,
295 : const casacore::Float& minbeamfrac=0.0,
296 : const casacore::Int growiterations=0,
297 : const casacore::Bool dogrowprune=true,
298 : const casacore::Float& minpercentchange=0.0,
299 : const casacore::Bool verbose=false,
300 : const casacore::Bool fastnoise=false,
301 : const casacore::Bool isthresholdreached=false,
302 : casacore::Float pblimit=0.1);
303 :
304 :
305 : // depth-first-search algorithm for 2D
306 : void depthFirstSearch(casacore::Int x,
307 : casacore::Int y,
308 : casacore::Int cur_label,
309 : casacore::Array<casacore::Float>& inlatarr,
310 : casacore::Array<casacore::Float>& lablatarr);
311 :
312 : // non-recursive depth-first-search algorithm for 2D
313 : void depthFirstSearch2(casacore::Int x,
314 : casacore::Int y,
315 : casacore::Int cur_label,
316 : casacore::Array<casacore::Float>& inlatarr,
317 : casacore::Array<casacore::Float>& lablatarr);
318 :
319 : // returns a Vector of neighboring pixels in IPosition (4-direction connectivity)
320 : casacore::Vector<casacore::IPosition> defineNeighbors(casacore::IPosition& pos,
321 : casacore::Int nrow,
322 : casacore::Int ncol);
323 :
324 : // label connected regions using depth-first-search algorithm
325 : void labelRegions(casacore::Lattice<casacore::Float>& inlat, casacore::Lattice<casacore::Float>& lablat);
326 :
327 :
328 :
329 : // find sizes of bolbs (regions) found by labelRegions
330 : casacore::Vector<casacore::Float> findBlobSize(casacore::Lattice<casacore::Float>& lablat);
331 :
332 : // check if mask image is empty (all zeros ) =True or not
333 : casacore::Bool isEmptyMask(casacore::ImageInterface<casacore::Float>& maskiamge);
334 : casacore::Int getTotalPixels(casacore::ImageInterface<casacore::Float>& maskiamge);
335 :
336 : // for warning messages for empy initial mask in automask
337 : //void noMaskCheck(casacore::ImageInterface<casacore::Float>& mask, casacore::Vector<casacore::String>& thresholdType);
338 : void noMaskCheck(casacore::ImageInterface<casacore::Float>& mask, casacore::Matrix<casacore::String>& thresholdType);
339 :
340 : // determining skip channels for the mask changed less than the specfied percentage
341 : void skipChannels(const casacore::Float& fracChnage,
342 : casacore::ImageInterface<casacore::Float>& prevmask,
343 : casacore::ImageInterface<casacore::Float>& curmask,
344 : const casacore::Matrix<casacore::String>& threshtype,
345 : const casacore::Bool isthresholdreached,
346 : casacore::Vector<casacore::Bool>& chanFlag,
347 : casacore::Vector<casacore::Bool>& zeroChanMask);
348 :
349 : // check if input image is a mask image with 0 or a value (if normalize=true, 1)
350 : //casacore::Bool checkMaskImage(casacore::ImageInterface<casacore::Float>& maskiamge, casacore::Bool normalize=true);
351 :
352 : // print per-channel automask summary
353 : void printAutomaskSummary(const casacore::Array<casacore::Double>& rmss,
354 : const casacore::Array<casacore::Double>& maxs,
355 : const casacore::Array<casacore::Double>& mins,
356 : const casacore::Array<casacore::Double>& mdns,
357 : //const casacore::Vector<casacore::Float>& maskthreshold,
358 : const casacore::Matrix<casacore::Float>& maskthreshold,
359 : //const casacore::Vector<casacore::String>& masktype,
360 : const casacore::Matrix<casacore::String>& masktype,
361 : const casacore::Vector<casacore::Bool>& chanflag,
362 : //const casacore::Vector<casacore::Bool>& zerochanmask,
363 : const casacore::Matrix<casacore::Bool>& zerochanmask,
364 : //const casacore::Vector<casacore::uInt>& nreg,
365 : const casacore::Matrix<casacore::uInt>& nreg,
366 : //const casacore::Vector<casacore::uInt>& npruned,
367 : const casacore::Matrix<casacore::uInt>& npruned,
368 : //const casacore::Vector<casacore::uInt>& ngrowreg,
369 : const casacore::Matrix<casacore::uInt>& ngrowreg,
370 : //const casacore::Vector<casacore::uInt>& ngrowpruned,
371 : const casacore::Matrix<casacore::uInt>& ngrowpruned,
372 : //const casacore::Vector<casacore::Float>& negmaskpixs,
373 : const casacore::Matrix<casacore::Float>& negmaskpixs,
374 : const casacore::Record& miscsummaryinfo);
375 :
376 :
377 : //
378 : casacore::Bool compareSpectralCoordinate(const casacore::ImageInterface<casacore::Float>& inImage,
379 : const casacore::ImageInterface<casacore::Float>& outImage);
380 : static casacore::Bool cloneImShape(const casacore::ImageInterface<casacore::Float>& inImage, const casacore::String& outImageName);
381 : // max MB of memory to use in TempImage
382 0 : static inline casacore::Double memoryToUse() {return -1.0;};
383 :
384 : // Calculate statistics on a residual image with additional region and LEL mask specifications
385 : // using classical method
386 : static casacore::Record calcImageStatistics(casacore::ImageInterface<casacore::Float>& res,
387 : casacore::String& lelmask,
388 : casacore::Record* regionPtr,
389 : const casacore::Bool robust);
390 :
391 : // Calcuate statistics on a residual image using robust methods to estimate RMS noise.
392 : // Basing on presence of a mask the following logic is used.
393 : // a. If there is no existing clean mask, calculate statistics using the chauvenet algorithm with maxiter=5 and zscore=-1.
394 : // b. If there is an existing clean mask, calculate the classic statistics with robust=True in the region outside the clean mask
395 : // and inside the primary beam mask.
396 : static casacore::Record calcRobustImageStatisticsOld(casacore::ImageInterface<casacore::Float>& res,
397 : casacore::ImageInterface<casacore::Float>& prevmask,
398 : casacore::LatticeExpr<casacore::Bool>& pbmask,
399 : casacore::String& lelmask,
400 : casacore::Record* regionPtr,
401 : const casacore::Bool robust,
402 : casacore::Vector<casacore::Bool>& chanflag);
403 :
404 : static casacore::Record calcRobustImageStatistics(casacore::ImageInterface<casacore::Float>& res,
405 : casacore::ImageInterface<casacore::Float>& prevmask,
406 : casacore::LatticeExpr<casacore::Bool>& pbmask,
407 : casacore::String& lelmask,
408 : casacore::Record* regionPtr,
409 : const casacore::Bool robust,
410 : casacore::Vector<casacore::Bool>& chanflag);
411 :
412 : // Store pbmask level (a.k.a pblimit for mask)
413 : void setPBMaskLevel(const casacore::Float pbmasklevel);
414 : casacore::Float getPBMaskLevel();
415 :
416 : private:
417 : double itsRms;
418 : double itsMax;
419 : float itsSidelobeLevel;
420 : float itsPBMaskLevel;
421 : int itsNAME_MAX;
422 : bool itsTooLongForFname;
423 : };
424 :
425 :
426 :
427 : } //# NAMESPACE CASA - END
428 :
429 :
430 : #endif
431 :
432 :
|
