Line data Source code
1 : //# ImageMoments.h: generate moments from images
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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: ImageMoments.h 20299 2008-04-03 05:56:44Z gervandiepen $
27 :
28 : #ifndef IMAGES_IMAGEMOMENTS_H
29 : #define IMAGES_IMAGEMOMENTS_H
30 :
31 : #include <imageanalysis/ImageAnalysis/MomentsBase.h>
32 :
33 : #include <imageanalysis/ImageTypedefs.h>
34 :
35 : namespace casacore{
36 :
37 : template <class T> class MaskedLattice;
38 : }
39 :
40 : namespace casa {
41 :
42 : class ImageMomentsProgressMonitor;
43 :
44 : // <summary>
45 : // This class generates moments from an image.
46 : // </summary>
47 :
48 : // <use visibility=export>
49 :
50 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
51 : // </reviewed>
52 :
53 : // <prerequisite>
54 : // <li> <linkto class="casacore::ImageInterface">casacore::ImageInterface</linkto>
55 : // <li> <linkto class="MomentsBase">MomentsBase</linkto>
56 : // <li> <linkto class="casacore::LatticeApply">casacore::LatticeApply</linkto>
57 : // <li> <linkto class="MomentCalcBase">MomentCalcBase</linkto>
58 : // </prerequisite>
59 :
60 : // <etymology>
61 : // This class computes moments from images
62 : // </etymology>
63 :
64 : // <synopsis>
65 : // The primary goal of this class is to help spectral-line astronomers analyze
66 : // their multi-dimensional images by generating moments of a specified axis.
67 : // The word "moment" is used loosely here. It refers to collapsing an axis
68 : // to one pixel and putting the value of that pixel (for all of the other
69 : // non-collapsed axes) to something computed from the data values along
70 : // the moment axis. For example, take an RA-DEC-Velocity cube, collapse
71 : // the velocity axis by computing the mean intensity at each RA-DEC
72 : // pixel. This class offers many different moments and a variety of
73 : // interactive and automatic ways to compute them.
74 : //
75 : // This class only accepts images of type <src>casacore::Float</src> and <src>casacore::Double</src>.
76 : // This restriction is because of the plotting capabilities which are a
77 : // bit awkward for other types.
78 : //
79 : // This class makes a distinction between a "moment" and a "method". This
80 : // boundary is a little blurred, but it claims to refer to the distinction
81 : // of what you are computing, as to how the pixels that were included in the
82 : // computation were selected. For example, a "moment" would be the average value
83 : // of some pixels. A method for selecting those pixels would be a simple range
84 : // specifying a range for which pixels should be included.
85 : //
86 : // The default state of this class is to do nothing. If you specify an image via
87 : // the <src>setImage</src> function then invoking the <src>createMoments</src>
88 : // function will cause it to compute the integrated itensity along the
89 : // spectral axis of the image (if it can find one). You can change the
90 : // computational state of the class from this basic form via the remaining
91 : // <src>set</src> functions. You can call any number of these functions to
92 : // suit your needs.
93 : //
94 : // Because there are a wide variety of methods available, if you specify an
95 : // invalid combination, a table showing the available methods is listed. It
96 : // is reproduced below for convenience.
97 : //
98 : // The basic method is to just compute moments directly from the pixel intensities.
99 : // This can be modified by applying pixel intensity inclusion or exclusion ranges.
100 : // You can then also smooth the image and compute a mask based on the inclusion or
101 : // exclusion ranges applied to the smoothed image. This mask is then applied to
102 : // the unsmoothed data for moment computation.
103 : //
104 : // The window method does no pixel intensity range selection. Instead a spectral
105 : // window is found (hopefully surrounding the spectral line feature) and only the
106 : // pixels in that window are used for computation. This window can be found from the
107 : // smoothed or unsmoothed data. The moments are always computed from the unsmoothed
108 : // data. For each spectrum, the window can be found interactively or automatically.
109 : // There are two interactive methods. Either you just mark the window with the
110 : // cursor, or you interactively fit a Gaussian to the profile and the +/- 3-sigma
111 : // window is returned. There are two automatic methods. Either Bosma's converging
112 : // mean algorithm is used, or an automatically fit Gaussian +/- 3-sigma window
113 : // is returned.
114 : //
115 : // The fitting method fits Gaussians to spectral features either automatically
116 : // or interactively. The moments are then computed from the Gaussian fits
117 : // (not the data themselves).
118 : //
119 : // When an output image is created, it will have N-1 axes, where the input image
120 : // has N axes. In the output image, the physical axis corresponding to the moment
121 : // axis will have been removed, but the coordinate information will be retained
122 : // for future coordinate transformations. For example, if you have a RA-DEC-VELOCITY
123 : // image and you collapsed axis 2 (the DEC axis) the output images would be
124 : // RA-VELOCITY with the coordinate information retained for the DEC axis so that
125 : // the coupled nature of RA/DEC coordinates is preserved.
126 : //
127 : // Output images are created with an all true (good) mask. If, for a given
128 : // pixel, the moment calculation fails, then the mask is set to F.
129 : //
130 : // When making plots, the order in which the spectra are displayed is determined
131 : // by the tiling sequence of the image (for optimum speed of access).
132 : //
133 : //
134 : // <srcblock>
135 : // Allowed Methods
136 : // ---------------
137 : //
138 : // casacore::Smooth Window Fit in/exclude Interactive
139 : // -----------------------------------------------------
140 : // N N N N N
141 : // Y/N N N Y N
142 : //
143 : // Y/N Y N N Y
144 : // Y/N Y N N N
145 : // Y/N Y Y N Y/N
146 : //
147 : // N N Y N Y/N
148 : // ----------------------------------------------------
149 : // </srcblock>
150 : //
151 : // <note role=caution>
152 : // Note that the <src>MEDIAN_COORDINATE</src> moment is not very robust.
153 : // It is very useful for generating quickly a velocity field in a way
154 : // that is not sensitive to noise. However, it will only give sensible
155 : // results under certain conditions. It treats the spectrum as a
156 : // probability distribution, generates the cumulative distribution for
157 : // the selected pixels (via an <src>include</src> or <src>exclude</src>
158 : // pixel range, and finds the (interpolated) coordinate coresponding to
159 : // the 50% value. The generation of the cumulative distribution and the
160 : // finding of the 50% level really only makes sense if the cumulative
161 : // distribution is monotonically increasing. This essentially means only
162 : // selecting pixels which are positive or negative. For this reason, this
163 : // moment type is *only* supported with the basic method (i.e. no smoothing,
164 : // no windowing, no fitting) with a pixel selection range that is either
165 : // all positive, or all negative.
166 : // </note>
167 : //
168 : // <note role=tip>
169 : // If you ignore return error statuses from the functions that set the
170 : // state of the class, the internal status of the class is set to bad.
171 : // This means it will just keep on returning error conditions until you
172 : // explicitly recover the situation. A message describing the last
173 : // error condition can be recovered with function errorMessage.
174 : // </note>
175 : // </synopsis>
176 :
177 : // <example>
178 : // <srcBlock>
179 : //// Set state function argument values
180 : //
181 : // casacore::Vector<casacore::Int> moments(2);
182 : // moments(0) = ImageMoments<casacore::Float>::AVERAGE;
183 : // moments(1) = ImageMoments<casacore::Float>::WEIGHTED_MEAN_COORDINATE;
184 : // casacore::Vector<int> methods(2);
185 : // methods(0) = ImageMoments<casacore::Float>::WINDOW;
186 : // methods(1) = ImageMoments<casacore::Float>::INTERACTIVE;
187 : // casacore::Vector<casacore::Int> nxy(2);
188 : // nxy(0) = 3;
189 : // nxy(1) = 3;
190 : //
191 : //// Open paged image
192 : //
193 : // casacore::PagedImage<casacore::Float> inImage(inName);
194 : //
195 : //// Construct moment helper object
196 : //
197 : // casacore::LogOrigin or("myClass", "myFunction(...)", WHERE);
198 : // casacore::LogIO os(or);
199 : // ImageMoments<casacore::Float> moment(casacore::SubImage<casacore::Float>(inName), os);
200 : //
201 : //// Specify state via control functions
202 : //
203 : // if (!moment.setMoments(moments)) return 1;
204 : // if (!moment.setWinFitMethod(methods)) return 1;
205 : // if (!moment.setMomentAxis(3)) return 1;
206 : // if (!moment.setPlotting("/xs", nxy)) return 1;
207 : //
208 : //// Create the moments
209 : //
210 : // if (!moment.createMoments()) return 1;
211 : // </srcBlock>
212 : // In this example, we generate two moments (average intensity and intensity
213 : // weighted mean coordinate -- usually the velocity field) of axis 3 by the
214 : // interactive window method. The interactive plotting is done on the
215 : // device called <src>/xs</src> (no longer supported). We put 9 subplots on each page. The output
216 : // file names are constructed by the class from the input file name plus some
217 : // internally generated suffixes.
218 : // </example>
219 :
220 : // <motivation>
221 : // This is a fundamental and traditional piece of spectral-line image analysis.
222 : // </motivation>
223 :
224 : // <todo asof="1996/11/26">
225 : // <li> more control over histogram of image noise at start (pixel
226 : // range and number of bins)
227 : // <li> better algorithm for seeing if spectrum is pure noise
228 : // <li> Make this class extensible so users could add their own method.
229 : // </todo>
230 :
231 :
232 : template <class T> class ImageMoments : public MomentsBase<T> {
233 : public:
234 :
235 : // Note that if I don't put MomentCalcBase as a forward declaration
236 : // and use instead "friend class MomentCalcBase<T>" The gnu compiler
237 : // fails with a typedef problem. But I can't solve it with say
238 : // <src>typedef MomentCalcBase<T> gpp_type;</src> because you can only do a
239 : // typedef with an actual type, not a <tt>T</tt> !
240 : friend class MomentCalcBase<T>;
241 :
242 : ImageMoments() = delete;
243 :
244 : // Constructor takes an image and a <src>casacore::LogIO</src> object for logging purposes.
245 : // You specify whether output images are automatically overwritten if pre-existing,
246 : // or whether an intercative choice dialog widget appears (overWriteOutput=F)
247 : // You may also determine whether a progress meter is displayed or not.
248 : ImageMoments (
249 : const casacore::ImageInterface<T>& image,
250 : casacore::LogIO &os,
251 : casacore::Bool overWriteOutput=false,
252 : casacore::Bool showProgress=true
253 : );
254 :
255 : ImageMoments(const ImageMoments<T> &other) = delete;
256 :
257 : // Destructor
258 : ~ImageMoments();
259 :
260 : ImageMoments<T> &operator=(const ImageMoments<T> &other) = delete;
261 :
262 : // Set the moment axis (0 relative). A return value of <src>false</src>
263 : // indicates that the axis was not contained in the image. If you don't
264 : // call this function, the default state of the class is to set the
265 : // moment axis to the spectral axis if it can find one. Otherwise
266 : // an error will result.
267 : casacore::Bool setMomentAxis (casacore::Int momentAxis);
268 :
269 : // This function invokes smoothing of the input image. Give <src>casacore::Int</src>
270 : // arrays for the axes (0 relative) to be smoothed and the smoothing kernel
271 : // types (use the <src>enum KernelTypes</src>) for each axis. Give a
272 : // <src>casacore::Double</src> array for the widths (full width for BOXCAR and full
273 : // width at half maximum for GAUSSIAN) in pixels of the smoothing kernels for
274 : // each axis. For HANNING smoothing, you always get the quarter-half-quarter
275 : // kernel (no matter what you might ask for). A return value of <src>false</src>
276 : // indicates that you have given an inconsistent or invalid set of smoothing
277 : // parameters. If you don't call this function the default state of the
278 : // class is to do no smoothing. The kernel types are specified with
279 : // the casacore::VectorKernel::KernelTypes enum
280 : casacore::Bool setSmoothMethod(
281 : const casacore::Vector<casacore::Int>& smoothAxes,
282 : const casacore::Vector<casacore::Int>& kernelTypes,
283 : const casacore::Vector<casacore::Quantum<casacore::Double> >& kernelWidths
284 : );
285 :
286 : casacore::Bool setSmoothMethod(
287 : const casacore::Vector<casacore::Int>& smoothAxes,
288 : const casacore::Vector<casacore::Int>& kernelTypes,
289 : const casacore::Vector<casacore::Double>& kernelWidths
290 : );
291 :
292 : // This is the function that does all the computational work. It should be called
293 : // after the <src>set</src> functions.
294 : // If the axis being collapsed comes from a coordinate with one axis only,
295 : // the axis and its coordinate are physically removed from the output image. Otherwise,
296 : // if <src>removeAxes=true</src> then the output axis is logically removed from the
297 : // the output CoordinateSystem. If <src>removeAxes=false</src> then the axis
298 : // is retained with shape=1 and with its original coordinate information (which
299 : // is probably meaningless).
300 : //
301 : // The output vector will hold PagedImages or TempImages (doTemp=true).
302 : // If doTemp is true, the outFileName is not used.
303 : //
304 : // If you create PagedImages, outFileName is the root name for
305 : // the output files. Suffixes will be made up internally to append
306 : // to this root. If you only ask for one moment,
307 : // this will be the actual name of the output file. If you don't set this
308 : // variable, the default state of the class is to set the output name root to
309 : // the name of the input file.
310 : std::vector<std::shared_ptr<casacore::MaskedLattice<T> > > createMoments(
311 : casacore::Bool doTemp, const casacore::String& outFileName,
312 : casacore::Bool removeAxes=true
313 : );
314 :
315 : // Set a new image. A return value of <src>false</src> indicates the
316 : // image had an invalid type (this class only accepts casacore::Float or casacore::Double images).
317 : casacore::Bool setNewImage (const casacore::ImageInterface<T>& image);
318 :
319 : // Get CoordinateSystem
320 0 : const casacore::CoordinateSystem& coordinates() {return _image->coordinates();};
321 :
322 : // Get shape
323 0 : casacore::IPosition getShape() const { return _image->shape(); }
324 :
325 : //Set an ImageMomentsProgressMonitor interested in getting updates on the
326 : //progress of the collapse process.
327 : void setProgressMonitor(ImageMomentsProgressMonitor* progressMonitor);
328 :
329 : private:
330 :
331 : SPCIIT _image = SPCIIT(nullptr);
332 : ImageMomentsProgressMonitor* _progressMonitor = nullptr;
333 :
334 : // casacore::Smooth an image
335 : SPIIT _smoothImage();
336 :
337 : // Determine the noise by fitting a Gaussian to a histogram
338 : // of the entire image above the 25% levels. If a plotting
339 : // device is set, the user can interact with this process.
340 : void _whatIsTheNoise (T& noise, const casacore::ImageInterface<T>& image);
341 :
342 : protected:
343 : using MomentsBase<T>::os_p;
344 : using MomentsBase<T>::showProgress_p;
345 : using MomentsBase<T>::momentAxisDefault_p;
346 : using MomentsBase<T>::peakSNR_p;
347 : using MomentsBase<T>::stdDeviation_p;
348 : using MomentsBase<T>::yMin_p;
349 : using MomentsBase<T>::yMax_p;
350 : using MomentsBase<T>::out_p;
351 : using MomentsBase<T>::smoothOut_p;
352 : using MomentsBase<T>::goodParameterStatus_p;
353 : using MomentsBase<T>::doWindow_p;
354 : using MomentsBase<T>::doFit_p;
355 : using MomentsBase<T>::doSmooth_p;
356 : using MomentsBase<T>::noInclude_p;
357 : using MomentsBase<T>::noExclude_p;
358 : using MomentsBase<T>::fixedYLimits_p;
359 : using MomentsBase<T>::momentAxis_p;
360 : using MomentsBase<T>::worldMomentAxis_p;
361 : using MomentsBase<T>::kernelTypes_p;
362 : using MomentsBase<T>::kernelWidths_p;
363 : using MomentsBase<T>::moments_p;
364 : using MomentsBase<T>::selectRange_p;
365 : using MomentsBase<T>::smoothAxes_p;
366 : using MomentsBase<T>::overWriteOutput_p;
367 : using MomentsBase<T>::error_p;
368 : using MomentsBase<T>::convertToVelocity_p;
369 : using MomentsBase<T>::velocityType_p;
370 : using MomentsBase<T>::_checkMethod;
371 : };
372 :
373 : }
374 :
375 : #ifndef CASACORE_NO_AUTO_TEMPLATES
376 : #include <imageanalysis/ImageAnalysis/ImageMoments.tcc>
377 : #endif
378 : #endif
|
