Line data Source code
1 : //# MomentCalcBase.h:
2 : //# Copyright (C) 1997,1999,2000,2001,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 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: MomentCalculator.h 20299 2008-04-03 05:56:44Z gervandiepen $
27 :
28 : #ifndef IMAGEANALYSIS_MOMENTCALCBASE_H
29 : #define IMAGEANALYSIS_MOMENTCALCBASE_H
30 :
31 : namespace casa {
32 :
33 : //# Forward Declarations
34 : template <class T> class MomentsBase;
35 : // <summary>
36 : // Abstract base class for moment calculator classes
37 : // </summary>
38 : // <use visibility=export>
39 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
40 : // </reviewed>
41 : //
42 : // <prerequisite>
43 : // <li> <linkto class="MomentsBase">MomentsBase</linkto>
44 : // <li> <linkto class="ImageMoments">ImageMoments</linkto>
45 : // <li> <linkto class="casacore::LatticeApply">casacore::LatticeApply</linkto>
46 : // <li> <linkto class="casacore::LineCollapser">casacore::LineCollapser</linkto>
47 : // </prerequisite>
48 : //
49 : // <synopsis>
50 : // This class, its concrete derived classes, and the classes casacore::LineCollapser,
51 : // ImageMoments, MSMoments, and casacore::LatticeApply are connected as follows. casacore::LatticeApply offers
52 : // functions so that the application programmer does not need to worry about how
53 : // to optimally iterate through a casacore::Lattice; it deals with tiling and to a
54 : // lesser extent memory. casacore::LatticeApply functions are used by offering a class
55 : // object to them that has a member function with a name and signature
56 : // specified by an abstract base class that casacore::LatticeApply uses and the
57 : // offered class inherits from. Specifically, in this case, MomentCalcBase
58 : // inherits from casacore::LineCollapser and casacore::LatticeApply uses objects and methods of this
59 : // class (but does not inherit from it). This defines the functions
60 : // <src>collapse</src> and <src>multiProcess</src> which operate on a vector
61 : // extracted from a Lattice. The former returns one number, the latter a vector
62 : // of numbers from that profile. MomentCalcBase is a base class for
63 : // for moment calculation and the <src>multiProcess</src>
64 : // functions are used to compute moments (e.g., mean, sum, sum squared,
65 : // intensity weighted velocity etc).
66 : //
67 : // It is actually the concrete classes derived from MomentCalcBase (call them,
68 : // as a group, the MomentCalculator classes) that implement the <src>multiProcess</src>
69 : // functions. These derived classes allow different
70 : // algorithms to be written with which moments of the vector can be computed.
71 : //
72 : // Now, so far, we have a casacore::LatticeApply function which iterates through Lattices,
73 : // extracts vectors, and offers them up to functions implemented in the derived
74 : // MomentCalculator classes to compute the moments. As well as that, we need some
75 : // class to actually construct the MomentCalculator classes and to feed them to
76 : // LatticeApply. This is the role of the ImageMoments or MSMoments classes.
77 : // They are a high level
78 : // class which takes control information from users specifying which moments they
79 : // would like to calculate and how. They also provide the ancilliary masking lattice to
80 : // the MomentCalculator constructors. The actual computational work is done by the
81 : // MomentCalculator classes. So MomentsBase, MomentCalcBase and their derived
82 : // MomentCalculator classes are really one unit; none of them are useful without
83 : // the others. The separation of functionality is caused by having the
84 : // casacore::LatticeApply class that knows all about optimally iterating through Lattices.
85 : //
86 : // The coupling between these classes is done partly by the "friendship". MomentsBase and
87 : // its inheritances
88 : // grant friendship to MomentCalcBase so that the latter has access to the private data and
89 : // private functions of the formers. MomentCalcBase then operates as an interface between
90 : // its derived MomentCalculator classes and ImageMoments or MSMoments. It retrieves private data
91 : // from these classes, and also activates private functions in these classes, on behalf
92 : // of the MomentCalculator classes. The rest of the coupling is done via the constructors
93 : // of the derived MomentCalculator classes.
94 : //
95 : // Finally, MomentCalcBase also has a number of protected functions that are common to its
96 : // derived classes (e.g. fitting, accumulating sums etc). It also has protected
97 : // data that is common to all the MomentCalculator classes. This protected data is accessed
98 : // directly by name rather than with interface functions as there is too much of it. Of
99 : // course, since MomentCalcBase is an abstract base class, it is up to the MomentCalculator
100 : // classes to give the MomentCalcBase protected data objects values.
101 : //
102 : // For discussion about different moments and algorithms to compute them see the
103 : // discussion in <linkto class="MomentsBase">MomentsBase</linkto>,
104 : // <linkto class="ImageMoments">ImageMoments</linkto>,
105 : // <linkto class="MSMoments">MSMoments</linkto> and also in
106 : // the derived classes documentation.
107 : // </synopsis>
108 : //
109 : // <example>
110 : // Since MomentCalcBase is an abstract class, we defer code examples to
111 : // the derived classes.
112 : // </example>
113 : //
114 : // <motivation>
115 : // We were desirous of writing functions to optimally iterate through Lattices
116 : // so that the application programmer did not have to know anything about tiling
117 : // or memory if possible. These are the casacore::LatticeApply functions. To incorporate
118 : // MomentsBase and its inheritances into this scheme required some of it to be shifted into
119 : // MomentCalcBase and its derived classes.
120 : // </motivation>
121 : //
122 : // <note role=tip>
123 : // Note that there are is assignment operator or copy constructor.
124 : // Do not use the ones the system would generate either.
125 : // </note>
126 : //
127 : // <todo asof="yyyy/mm/dd">
128 : // <li> Derive more classes !
129 : // </todo>
130 :
131 :
132 : template <class T> class MomentCalcBase : public casacore::LineCollapser<T,T> {
133 : public:
134 :
135 : using AccumType = typename casacore::NumericTraits<T>::PrecisionType;
136 : using DataIterator = typename casacore::Vector<T>::const_iterator;
137 : using MaskIterator = casacore::Vector<casacore::Bool>::const_iterator;
138 :
139 : virtual ~MomentCalcBase();
140 :
141 : // Returns the number of failed fits if doing fitting
142 0 : virtual inline casacore::uInt nFailedFits() const { return nFailed_p; }
143 :
144 : protected:
145 :
146 : // A number of private data members are kept here in the base class
147 : // as they are common to the derived classes. Since this class
148 : // is abstract, they have to be filled by the derived classes.
149 :
150 : // CoordinateSystem
151 : casacore::CoordinateSystem cSys_p;
152 :
153 : // This vector is a container for all the possible moments that
154 : // can be calculated. They are in the order given by the MomentsBase
155 : // enum MomentTypes
156 : casacore::Vector<T> calcMoments_p;
157 : casacore::Vector<casacore::Bool> calcMomentsMask_p;
158 :
159 : // This vector tells us which elements of the calcMoments_p vector
160 : // we wish to select
161 : casacore::Vector<casacore::Int> selectMoments_p;
162 :
163 : // Although the general philosophy of these classes is to compute
164 : // all the posisble moments and then select the ones we want,
165 : // some of them are too expensive to calculate unless they are
166 : // really wanted. These are the median moments and those that
167 : // require a second pass. These control Bools tell us whether
168 : // we really want to compute the expensive ones.
169 : casacore::Bool doMedianI_p, doMedianV_p, doAbsDev_p;
170 :
171 : // These vectors are used to transform coordinates between pixel and world
172 : casacore::Vector<casacore::Double> pixelIn_p, worldOut_p;
173 :
174 : // All computations involving casacore::Coordinate conversions are relatively expensive
175 : // These Bools signifies whether we need coordinate calculations or not for
176 : // the full profile, and for some occaisional calculations
177 : casacore::Bool doCoordProfile_p, doCoordRandom_p;
178 :
179 : // This vector houses the world coordinate values for the profile if it
180 : // was from a separable axis. This means this vector can be pre computed
181 : // just once, instead of working out the coordinates for each profile
182 : // (expensive). It should only be filled if doCoordCalc_p is true
183 : casacore::Vector<casacore::Double> sepWorldCoord_p;
184 :
185 : // This vector is used to hold the abscissa values
186 : casacore::Vector<T> abcissa_p;
187 :
188 : // This string tells us the name of the moment axis (VELO or FREQ etc)
189 : casacore::String momAxisType_p;
190 :
191 : // This is the number of Gaussian fits that failed.
192 : casacore::uInt nFailed_p;
193 :
194 : // This scale factor is the increment along the moment axis
195 : // applied so that units for the Integrated moment are like
196 : // Jy/beam.km/s (or whatever is needed for the moment axis units)
197 : // For non-linear velocities (e.g. optical) this is approximate
198 : // only and is computed at the reference pixel
199 : casacore::Double integratedScaleFactor_p;
200 :
201 : // Accumulate statistical sums from a vector
202 0 : inline void accumSums(
203 : typename casacore::NumericTraits<T>::PrecisionType& s0,
204 : typename casacore::NumericTraits<T>::PrecisionType& s0Sq,
205 : typename casacore::NumericTraits<T>::PrecisionType& s1,
206 : typename casacore::NumericTraits<T>::PrecisionType& s2,
207 : casacore::Int& iMin, casacore::Int& iMax, T& dMin, T& dMax,
208 : casacore::Int i, T datum, casacore::Double coord
209 : ) const {
210 : // Accumulate statistical sums from this datum
211 : //
212 : // casacore::Input:
213 : // i Index
214 : // datum Pixel value
215 : // coord casacore::Coordinate value on moment axis
216 : // casacore::Input/output:
217 : // iMin,max index of dMin and dMax
218 : // dMin,dMax minimum and maximum value
219 : // Output:
220 : // s0 sum (I)
221 : // s0Sq sum (I*I)
222 : // s1 sum (I*v)
223 : // s2 sum (I*v*v)
224 0 : typename casacore::NumericTraits<T>::PrecisionType dDatum = datum;
225 0 : s0 += dDatum;
226 0 : s0Sq += dDatum*dDatum;
227 0 : s1 += dDatum*coord;
228 0 : s2 += dDatum*coord*coord;
229 0 : if (datum < dMin) {
230 0 : iMin = i;
231 0 : dMin = datum;
232 : }
233 0 : if (datum > dMax) {
234 0 : iMax = i;
235 0 : dMax = datum;
236 : }
237 0 : }
238 :
239 : // Determine if the spectrum is pure noise
240 : casacore::uInt allNoise(T& dMean,
241 : const casacore::Vector<T>& data,
242 : const casacore::Vector<casacore::Bool>& mask,
243 : T peakSNR,
244 : T stdDeviation
245 : ) const;
246 :
247 : // Check validity of constructor inputs
248 : void constructorCheck(
249 : casacore::Vector<T>& calcMoments,
250 : casacore::Vector<casacore::Bool>& calcMomentsMask,
251 : const casacore::Vector<casacore::Int>& selectMoments,
252 : casacore::uInt nLatticeOut
253 : ) const;
254 :
255 : // Find out from the selectMoments array whether we want
256 : // to compute the more expensive moments
257 : void costlyMoments(
258 : MomentsBase<T>& iMom, casacore::Bool& doMedianI,
259 : casacore::Bool& doMedianV, casacore::Bool& doAbsDev
260 : ) const;
261 :
262 : // Return the casacore::Bool saying whether we need to compute coordinates
263 : // or not for the requested moments
264 : void doCoordCalc(
265 : casacore::Bool& doCoordProfile,
266 : casacore::Bool& doCoordRandom,
267 : const MomentsBase<T>& iMom
268 : ) const;
269 :
270 : // Return the casacore::Bool from the ImageMoments or MSMoments object saying whether we
271 : // are going to fit Gaussians to the profiles or not.
272 : casacore::Bool doFit(const MomentsBase<T>& iMom) const;
273 :
274 : // Find the next masked or unmasked point in a vector
275 : casacore::Bool findNextDatum(
276 : casacore::uInt& iFound, const casacore::uInt& n,
277 : const casacore::Vector<casacore::Bool>& mask, const casacore::uInt& iStart,
278 : const casacore::Bool& findGood
279 : ) const;
280 :
281 : // Fit a Gaussian to x and y arrays given guesses for the gaussian parameters
282 : casacore::Bool fitGaussian(
283 : casacore::uInt& nFailed, T& peak, T& pos,
284 : T& width, T& level, const casacore::Vector<T>& x,
285 : const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask,
286 : T peakGuess, T posGuess, T widthGuess,
287 : T levelGuess
288 : ) const;
289 :
290 : // Automatically fit a Gaussian to a spectrum, including finding the
291 : // starting guesses.
292 : casacore::Bool getAutoGaussianFit(
293 : casacore::uInt& nFailed, casacore::Vector<T>& gaussPars,
294 : const casacore::Vector<T>& x, const casacore::Vector<T>& y,
295 : const casacore::Vector<casacore::Bool>& mask, T peakSNR,
296 : T stdDeviation
297 : ) const;
298 :
299 : // Automatically work out a guess for the Gaussian parameters
300 : // Returns false if all pixels masked.
301 : casacore::Bool getAutoGaussianGuess(
302 : T& peakGuess, T& posGuess,
303 : T& widthGuess, T& levelGuess,
304 : const casacore::Vector<T>& x, const casacore::Vector<T>& y,
305 : const casacore::Vector<casacore::Bool>& mask
306 : ) const;
307 :
308 : // Compute the world coordinate for the given moment axis pixel
309 0 : inline casacore::Double getMomentCoord(
310 : const MomentsBase<T>& iMom, casacore::Vector<casacore::Double>& pixelIn,
311 : casacore::Vector<casacore::Double>& worldOut, casacore::Double momentPixel,
312 : casacore::Bool asVelocity=false
313 : ) const {
314 : // Find the value of the world coordinate on the moment axis
315 : // for the given moment axis pixel value.
316 : //
317 : // Input
318 : // momentPixel is the index in the profile extracted from the data
319 : // casacore::Input/output
320 : // pixelIn Pixels to convert. Must all be filled in except for
321 : // pixelIn(momentPixel).
322 : // worldOut casacore::Vector to hold result
323 : //
324 : // Should really return a casacore::Fallible as I don't check and see
325 : // if the coordinate transformation fails or not
326 : //
327 : // Should really check the result is true, but for speed ...
328 0 : pixelIn[iMom.momentAxis_p] = momentPixel;
329 0 : cSys_p.toWorld(worldOut, pixelIn);
330 0 : if (asVelocity) {
331 : casacore::Double velocity;
332 0 : cSys_p.spectralCoordinate().frequencyToVelocity(
333 0 : velocity, worldOut(iMom.worldMomentAxis_p)
334 : );
335 0 : return velocity;
336 : }
337 0 : return worldOut(iMom.worldMomentAxis_p);
338 : }
339 :
340 : // Examine a mask and determine how many segments of unmasked points
341 : // it consists of.
342 : void lineSegments (
343 : casacore::uInt& nSeg, casacore::Vector<casacore::uInt>& start,
344 : casacore::Vector<casacore::uInt>& nPts, const casacore::Vector<casacore::Bool>& mask
345 : ) const;
346 :
347 : // Return the moment axis from the ImageMoments object
348 : casacore::Int& momentAxis(MomentsBase<T>& iMom) const;
349 :
350 : // Return the name of the moment/profile axis
351 : casacore::String momentAxisName(
352 : const casacore::CoordinateSystem&,
353 : const MomentsBase<T>& iMom
354 : ) const;
355 :
356 : // Return the peak SNR for determination of all noise spectra from
357 : // the ImageMoments or MSMoments object
358 : T& peakSNR(MomentsBase<T>& iMom) const;
359 :
360 : // Return the selected pixel intensity range from the ImageMoments or MSMoments
361 : // object and the Bools describing whether it is inclusion or exclusion
362 : void selectRange(
363 : casacore::Vector<T>& pixelRange,
364 : casacore::Bool& doInclude,
365 : casacore::Bool& doExlude,
366 : MomentsBase<T>& iMom
367 : ) const;
368 :
369 : // The MomentCalculators compute a vector of all possible moments.
370 : // This function returns a vector which selects the desired moments from that
371 : // "all moment" vector.
372 : casacore::Vector<casacore::Int> selectMoments(MomentsBase<T>& iMom) const;
373 :
374 : // Fill the ouput moments array
375 : void setCalcMoments(
376 : const MomentsBase<T>& iMom, casacore::Vector<T>& calcMoments,
377 : casacore::Vector<casacore::Bool>& calcMomentsMask, casacore::Vector<casacore::Double>& pixelIn,
378 : casacore::Vector<casacore::Double>& worldOut, casacore::Bool doCoord,
379 : casacore::Double integratedScaleFactor, T dMedian,
380 : T vMedian, casacore::Int nPts,
381 : typename casacore::NumericTraits<T>::PrecisionType s0,
382 : typename casacore::NumericTraits<T>::PrecisionType s1,
383 : typename casacore::NumericTraits<T>::PrecisionType s2,
384 : typename casacore::NumericTraits<T>::PrecisionType s0Sq,
385 : typename casacore::NumericTraits<T>::PrecisionType sumAbsDev,
386 : T dMin, T dMax, casacore::Int iMin, casacore::Int iMax
387 : ) const;
388 :
389 : // Fill a string with the position of the cursor
390 : void setPosLabel(casacore::String& title, const casacore::IPosition& pos) const;
391 :
392 : // Install casacore::CoordinateSystem and SpectralCoordinate
393 : // in protected data members
394 : void setCoordinateSystem(
395 : casacore::CoordinateSystem& cSys, MomentsBase<T>& iMom
396 : );
397 :
398 : // Set up separable moment axis coordinate vector and
399 : // conversion vectors if not separable
400 : void setUpCoords(
401 : const MomentsBase<T>& iMom, casacore::Vector<casacore::Double>& pixelIn,
402 : casacore::Vector<casacore::Double>& worldOut, casacore::Vector<casacore::Double>& sepWorldCoord,
403 : casacore::LogIO& os, casacore::Double& integratedScaleFactor,
404 : const casacore::CoordinateSystem& cSys, casacore::Bool doCoordProfile,
405 : casacore::Bool doCoordRandom
406 : ) const;
407 :
408 : // Return standard deviation of image from ImageMoments or MSMoments object
409 : T& stdDeviation(MomentsBase<T>& iMom) const;
410 :
411 : protected:
412 : // Check if #pixels is indeed 1.
413 : virtual void init (casacore::uInt nOutPixelsPerCollapse);
414 : };
415 :
416 : }
417 :
418 : #ifndef CASACORE_NO_AUTO_TEMPLATES
419 : #include <imageanalysis/ImageAnalysis/MomentCalcBase.tcc>
420 : #endif
421 : #endif
|
