Line data Source code
1 : //# MomentCalcBase.tcc:
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2004
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: MomentCalculator.tcc 19940 2007-02-27 05:35:22Z Malte.Marquarding $
27 : //
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/Arrays/Vector.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : #include <casacore/scimath/Fitting/NonLinearFitLM.h>
35 : #include <casacore/scimath/Functionals/Polynomial.h>
36 : #include <casacore/scimath/Functionals/CompoundFunction.h>
37 : #include <imageanalysis/ImageAnalysis/MomentsBase.h>
38 : #include <casacore/lattices/LatticeMath/LatticeStatsBase.h>
39 : #include <casacore/casa/BasicMath/Math.h>
40 : #include <casacore/casa/Logging/LogIO.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/casa/Exceptions/Error.h>
43 :
44 : namespace casa {
45 :
46 0 : template <class T> MomentCalcBase<T>::~MomentCalcBase() {}
47 :
48 0 : template <class T> void MomentCalcBase<T>::init(
49 : casacore::uInt nOutPixelsPerCollapse
50 : ) {
51 0 : AlwaysAssert (nOutPixelsPerCollapse == 1, casacore::AipsError);
52 0 : }
53 :
54 0 : template <class T> casacore::uInt MomentCalcBase<T>::allNoise (
55 : T& dMean, const casacore::Vector<T>& data, const casacore::Vector<casacore::Bool>& mask,
56 : const T peakSNR, const T stdDeviation
57 : ) const {
58 : // Try and work out whether this spectrum is all noise
59 : // or not. We don't bother with it if it is noise.
60 : // We compare the peak with sigma and a cutoff SNR
61 : // Returns 1 if all noise
62 : // Returns 2 if all masked
63 : // Returns 0 otherwise
64 0 : casacore::ClassicalStatistics<AccumType, DataIterator, MaskIterator> statsCalculator;
65 0 : statsCalculator.addData(data.begin(), mask.begin(), data.size());
66 0 : StatsData<AccumType> stats = statsCalculator.getStatistics();
67 0 : if (stats.npts == 0) {
68 0 : return 2;
69 : }
70 0 : T dMin = *stats.min;
71 0 : T dMax = *stats.max;
72 0 : dMean = stats.mean;
73 : // Assume we are continuum subtracted so outside of line mean=0
74 0 : const T rat = max(abs(dMin),abs(dMax)) / stdDeviation;
75 0 : casacore::uInt ret = rat < peakSNR ? 1 : 0;
76 0 : return ret;
77 0 : }
78 :
79 : template <class T>
80 0 : void MomentCalcBase<T>::constructorCheck(casacore::Vector<T>& calcMoments,
81 : casacore::Vector<casacore::Bool>& calcMomentsMask,
82 : const casacore::Vector<casacore::Int>& selectMoments,
83 : const casacore::uInt nLatticeOut) const
84 : {
85 : // Number of output lattices must equal the number of moments
86 : // the user asked to calculate
87 :
88 0 : AlwaysAssert(nLatticeOut == selectMoments.nelements(), casacore::AipsError);
89 :
90 : // Number of requested moments must be in allowed range
91 :
92 0 : auto nMaxMoments = MomentsBase<T>::NMOMENTS;
93 0 : AlwaysAssert(selectMoments.nelements() <= nMaxMoments, casacore::AipsError);
94 0 : AlwaysAssert(selectMoments.nelements() > 0, casacore::AipsError);
95 :
96 : // Resize the vector that will hold ALL possible moments
97 :
98 0 : calcMoments.resize(nMaxMoments);
99 0 : calcMomentsMask.resize(nMaxMoments);
100 0 : }
101 :
102 :
103 :
104 :
105 :
106 : template <class T>
107 0 : void MomentCalcBase<T>::costlyMoments(MomentsBase<T>& iMom,
108 : casacore::Bool& doMedianI,
109 : casacore::Bool& doMedianV,
110 : casacore::Bool& doAbsDev) const
111 : {
112 0 : doMedianI = false;
113 0 : doMedianV = false;
114 0 : doAbsDev = false;
115 : using IM = MomentsBase<casacore::Float>;
116 : //
117 0 : for (casacore::uInt i=0; i<iMom.moments_p.nelements(); i++) {
118 0 : if (iMom.moments_p(i) == IM::MEDIAN) doMedianI = true;
119 0 : if (iMom.moments_p(i) == IM::MEDIAN_COORDINATE) doMedianV = true;
120 0 : if (iMom.moments_p(i) == IM::ABS_MEAN_DEVIATION) doAbsDev = true;
121 : }
122 0 : }
123 :
124 : template <class T>
125 0 : Bool MomentCalcBase<T>::doFit(const MomentsBase<T>& iMom) const
126 : {
127 : // Get it from ImageMoments private data
128 :
129 0 : return iMom.doFit_p;
130 : }
131 :
132 : template <class T>
133 0 : void MomentCalcBase<T>::doCoordCalc(casacore::Bool& doCoordProfile,
134 : casacore::Bool& doCoordRandom,
135 : const MomentsBase<T>& iMom) const
136 : //
137 : // doCoordProfile - we need the coordinate for each pixel of the profile
138 : // doCoordRandom - we need the coordinate for occaisional use
139 : //
140 : {
141 : // Figure out if we need to compute the coordinate of each profile pixel index
142 : // for each profile. This is very expensive for non-separable axes.
143 :
144 0 : doCoordProfile = false;
145 0 : doCoordRandom = false;
146 : using IM = MomentsBase<casacore::Float>;
147 : //
148 0 : for (casacore::uInt i=0; i<iMom.moments_p.nelements(); i++) {
149 0 : if (iMom.moments_p(i) == IM::WEIGHTED_MEAN_COORDINATE ||
150 0 : iMom.moments_p(i) == IM::WEIGHTED_DISPERSION_COORDINATE) {
151 0 : doCoordProfile = true;
152 : }
153 0 : if (iMom.moments_p(i) == IM::MAXIMUM_COORDINATE ||
154 0 : iMom.moments_p(i) == IM::MINIMUM_COORDINATE ||
155 0 : iMom.moments_p(i) == IM::MEDIAN_COORDINATE) {
156 0 : doCoordRandom = true;
157 : }
158 : }
159 0 : }
160 :
161 : template <class T>
162 : Bool MomentCalcBase<T>::findNextDatum (casacore::uInt& iFound,
163 : const casacore::uInt& n,
164 : const casacore::Vector<casacore::Bool>& mask,
165 : const casacore::uInt& iStart,
166 : const casacore::Bool& findGood) const
167 : //
168 : // Find the next good (or bad) point in an array.
169 : // A good point in the array has a non-zero value.
170 : //
171 : // Inputs:
172 : // n Number of points in array
173 : // mask casacore::Vector containing counts.
174 : // iStart The index of the first point to consider
175 : // findGood If true look for next good point.
176 : // If false look for next bad point
177 : // Outputs:
178 : // iFound Index of found point
179 : // casacore::Bool false if didn't find another valid datum
180 : {
181 : for (casacore::uInt i=iStart; i<n; i++) {
182 : if ( (findGood && mask(i)) ||
183 : (!findGood && !mask(i)) ) {
184 : iFound = i;
185 : return true;
186 : }
187 : }
188 : return false;
189 : }
190 :
191 0 : template <class T> casacore::Bool MomentCalcBase<T>::fitGaussian(
192 : casacore::uInt& nFailed, T& peak, T& pos, T& width,
193 : T& level, const casacore::Vector<T>& x, const casacore::Vector<T>& y,
194 : const casacore::Vector<casacore::Bool>& mask, const T peakGuess,
195 : const T posGuess, const T widthGuess,
196 : const T levelGuess
197 : ) const {
198 : // Fit Gaussian pos * exp(-4ln2*(x-pos)**2/width**2)
199 : // width = fwhm
200 : // Returns false if fit fails or all masked
201 : // Select unmasked pixels
202 0 : casacore::uInt j = 0;
203 0 : auto nAll = y.size();
204 0 : casacore::Vector<T> xSel(nAll);
205 0 : casacore::Vector<T> ySel(nAll);
206 0 : for (casacore::uInt i=0; i<nAll; ++i) {
207 0 : if (mask[i]) {
208 0 : xSel[j] = x[i];
209 0 : ySel[j] = y[i];
210 0 : ++j;
211 : }
212 : }
213 0 : auto nPts = j;
214 0 : if (nPts == 0) {
215 0 : return false;
216 : }
217 0 : xSel.resize(nPts, true);
218 0 : ySel.resize(nPts, true);
219 : // Create fitter as gaussian + constant offset
220 0 : casacore::NonLinearFitLM<T> fitter;
221 0 : casacore::Gaussian1D<casacore::AutoDiff<T> > gauss;
222 0 : casacore::Polynomial<casacore::AutoDiff<T> > poly;
223 0 : casacore::CompoundFunction<casacore::AutoDiff<T> > func;
224 0 : func.addFunction(gauss);
225 0 : func.addFunction(poly);
226 0 : fitter.setFunction(func);
227 : // Initial guess
228 0 : casacore::Vector<T> v(4);
229 0 : v[0] = peakGuess;
230 0 : v[1] = posGuess;
231 0 : v[2] = widthGuess;
232 0 : v[3] = levelGuess;
233 0 : fitter.setParameterValues(v);
234 : // Set maximum number of iterations to 50. Default is 10
235 0 : fitter.setMaxIter(50);
236 : // Set converge criteria.
237 0 : fitter.setCriteria(0.001);
238 : // Perform fit on unmasked data
239 0 : casacore::Vector<T> resultSigma(nPts, 1);
240 0 : casacore::Vector<T> solution;
241 : try {
242 0 : solution = fitter.fit(xSel, ySel, resultSigma);
243 : }
244 0 : catch (const casacore::AipsError& x1) {
245 0 : ++nFailed;
246 0 : return false;
247 : }
248 : // Return values of fit
249 : // FIXME shouldn't these only be set if the fit converged?
250 0 : peak = solution[0];
251 0 : pos = solution[1];
252 0 : width = abs(solution[2]);
253 0 : level = solution[3];
254 : // Return status
255 0 : auto converged = fitter.converged();
256 0 : if (! converged) {
257 0 : ++nFailed;
258 : }
259 0 : return converged;
260 0 : }
261 :
262 : template <class T>
263 0 : Bool MomentCalcBase<T>::getAutoGaussianFit (casacore::uInt& nFailed,
264 : casacore::Vector<T>& gaussPars,
265 : const casacore::Vector<T>& x,
266 : const casacore::Vector<T>& y,
267 : const casacore::Vector<casacore::Bool>& mask,
268 : const T peakSNR,
269 : const T stdDeviation
270 : ) const
271 : //
272 : // Automatically fit a Gaussian and return the Gaussian parameters.
273 : // If a plotting device is active, we also plot the spectra and fits
274 : //
275 : // Inputs:
276 : // x,y casacore::Vector containing the data
277 : // mask true is good
278 : // plotter Plot spectrum and optionally the window
279 : // x,yLabel Labels
280 : // title
281 : // casacore::Input/output
282 : // nFailed Cumulative number of failed fits
283 : // Output:
284 : // gaussPars The gaussian parameters, peak, pos, fwhm
285 : // casacore::Bool If false then this spectrum has been rejected (all
286 : // masked, all noise, failed fit)
287 : //
288 : {
289 :
290 : // See if this spectrum is all noise. If so, forget it.
291 : // Return straight away if all masked
292 :
293 : T dMean;
294 0 : casacore::uInt iNoise = this->allNoise(dMean, y, mask, peakSNR, stdDeviation);
295 0 : if (iNoise == 2) return false;
296 :
297 0 : if (iNoise==1) {
298 0 : gaussPars = 0;
299 0 : return false;
300 : }
301 :
302 : // Work out guesses for Gaussian
303 :
304 : T peakGuess, posGuess, widthGuess, levelGuess;
305 : T pos, width, peak, level;
306 0 : if (!getAutoGaussianGuess(peakGuess, posGuess, widthGuess,
307 0 : levelGuess, x, y, mask)) return false;
308 0 : peakGuess = peakGuess - levelGuess;
309 :
310 :
311 : // Fit gaussian. Do it twice.
312 :
313 0 : if (!fitGaussian (nFailed, peak, pos, width, level, x, y, mask, peakGuess,
314 : posGuess, widthGuess, levelGuess)) {
315 0 : gaussPars = 0;
316 0 : return false;
317 : }
318 0 : gaussPars(0) = peak;
319 0 : gaussPars(1) = pos;
320 0 : gaussPars(2) = width;
321 0 : gaussPars(3) = level;
322 :
323 0 : return true;
324 : }
325 :
326 0 : template <class T> casacore::Bool MomentCalcBase<T>::getAutoGaussianGuess (
327 : T& peakGuess, T& posGuess, T& widthGuess, T& levelGuess,
328 : const casacore::Vector<T>& x, const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask
329 : ) const {
330 : // Make a wild stab in the dark as to what the Gaussian
331 : // parameters of this spectrum might be
332 :
333 0 : casacore::ClassicalStatistics<AccumType, DataIterator, MaskIterator> statsCalculator;
334 0 : statsCalculator.addData(y.begin(), mask.begin(), y.size());
335 0 : StatsData<AccumType> stats = statsCalculator.getStatistics();
336 0 : if (stats.npts == 0) {
337 : // all masked
338 0 : return false;
339 : }
340 : // Find peak and position of peak
341 0 : posGuess = x[stats.maxpos.second];
342 0 : peakGuess = *stats.max;
343 0 : levelGuess = stats.mean;
344 : // Nothing much is very robust. Assume the line is reasonably
345 : // sampled and set its width to a few pixels. Totally ridiculous.
346 0 : widthGuess = 5;
347 0 : return true;
348 0 : }
349 :
350 : template <class T>
351 : void MomentCalcBase<T>::lineSegments (casacore::uInt& nSeg,
352 : casacore::Vector<casacore::uInt>& start,
353 : casacore::Vector<casacore::uInt>& nPts,
354 : const casacore::Vector<casacore::Bool>& mask) const
355 : //
356 : // Examine an array and determine how many segments
357 : // of good points it consists of. A good point
358 : // occurs if the array value is greater than zero.
359 : //
360 : // Inputs:
361 : // mask The array mask. true is good.
362 : // Outputs:
363 : // nSeg Number of segments
364 : // start Indices of start of each segment
365 : // nPts Number of points in segment
366 : //
367 : {
368 : casacore::Bool finish = false;
369 : nSeg = 0;
370 : casacore::uInt iGood, iBad;
371 : const casacore::uInt n = mask.nelements();
372 : start.resize(n);
373 : nPts.resize(n);
374 :
375 : for (casacore::uInt i=0; !finish;) {
376 : if (!findNextDatum (iGood, n, mask, i, true)) {
377 : finish = true;
378 : } else {
379 : nSeg++;
380 : start(nSeg-1) = iGood;
381 :
382 : if (!findNextDatum (iBad, n, mask, iGood, false)) {
383 : nPts(nSeg-1) = n - start(nSeg-1);
384 : finish = true;
385 : } else {
386 : nPts(nSeg-1) = iBad - start(nSeg-1);
387 : i = iBad + 1;
388 : }
389 : }
390 : }
391 : start.resize(nSeg,true);
392 : nPts.resize(nSeg,true);
393 : }
394 :
395 : template <class T>
396 0 : Int& MomentCalcBase<T>::momentAxis(MomentsBase<T>& iMom) const
397 : {
398 0 : return iMom.momentAxis_p;
399 : }
400 :
401 : template <class T>
402 0 : String MomentCalcBase<T>::momentAxisName(const casacore::CoordinateSystem& cSys,
403 : const MomentsBase<T>& iMom) const
404 : {
405 : // Return the name of the moment/profile axis
406 :
407 0 : casacore::Int worldMomentAxis = cSys.pixelAxisToWorldAxis(iMom.momentAxis_p);
408 0 : return cSys.worldAxisNames()(worldMomentAxis);
409 : }
410 :
411 : template <class T>
412 0 : T& MomentCalcBase<T>::peakSNR(MomentsBase<T>& iMom) const
413 : {
414 : // Get it from ImageMoments private data
415 :
416 0 : return iMom.peakSNR_p;
417 : }
418 :
419 :
420 : template <class T>
421 0 : void MomentCalcBase<T>::selectRange(casacore::Vector<T>& pixelRange,
422 : casacore::Bool& doInclude,
423 : casacore::Bool& doExclude,
424 : MomentsBase<T>& iMom) const
425 : {
426 : // Get it from ImageMoments private data
427 :
428 0 : pixelRange = iMom.selectRange_p;
429 0 : doInclude = (!(iMom.noInclude_p));
430 0 : doExclude = (!(iMom.noExclude_p));
431 0 : }
432 :
433 :
434 : template <class T>
435 0 : Vector<casacore::Int> MomentCalcBase<T>::selectMoments(MomentsBase<T>& iMom) const
436 : //
437 : // Fill the moment selection vector according to what the user requests
438 : //
439 : {
440 : using IM = MomentsBase<casacore::Float>;
441 0 : casacore::Vector<casacore::Int> sel(IM::NMOMENTS);
442 :
443 0 : casacore::uInt j = 0;
444 0 : for (casacore::uInt i=0; i<iMom.moments_p.nelements(); i++) {
445 0 : if (iMom.moments_p(i) == IM::AVERAGE) {
446 0 : sel(j++) = IM::AVERAGE;
447 0 : } else if (iMom.moments_p(i) == IM::INTEGRATED) {
448 0 : sel(j++) = IM::INTEGRATED;
449 0 : } else if (iMom.moments_p(i) == IM::WEIGHTED_MEAN_COORDINATE) {
450 0 : sel(j++) = IM::WEIGHTED_MEAN_COORDINATE;
451 0 : } else if (iMom.moments_p(i) == IM::WEIGHTED_DISPERSION_COORDINATE) {
452 0 : sel(j++) = IM::WEIGHTED_DISPERSION_COORDINATE;
453 0 : } else if (iMom.moments_p(i) == IM::MEDIAN) {
454 0 : sel(j++) = IM::MEDIAN;
455 0 : } else if (iMom.moments_p(i) == IM::STANDARD_DEVIATION) {
456 0 : sel(j++) = IM::STANDARD_DEVIATION;
457 0 : } else if (iMom.moments_p(i) == IM::RMS) {
458 0 : sel(j++) = IM::RMS;
459 0 : } else if (iMom.moments_p(i) == IM::ABS_MEAN_DEVIATION) {
460 0 : sel(j++) = IM::ABS_MEAN_DEVIATION;
461 0 : } else if (iMom.moments_p(i) == IM::MAXIMUM) {
462 0 : sel(j++) = IM::MAXIMUM;
463 0 : } else if (iMom.moments_p(i) == IM::MAXIMUM_COORDINATE) {
464 0 : sel(j++) = IM::MAXIMUM_COORDINATE;
465 0 : } else if (iMom.moments_p(i) == IM::MINIMUM) {
466 0 : sel(j++) = IM::MINIMUM;
467 0 : } else if (iMom.moments_p(i) == IM::MINIMUM_COORDINATE) {
468 0 : sel(j++) = IM::MINIMUM_COORDINATE;
469 0 : } else if (iMom.moments_p(i) == IM::MEDIAN_COORDINATE) {
470 0 : sel(j++) = IM::MEDIAN_COORDINATE;
471 : }
472 : }
473 0 : sel.resize(j,true);
474 0 : return sel;
475 0 : }
476 :
477 :
478 : template <class T>
479 0 : void MomentCalcBase<T>::setPosLabel (casacore::String& title,
480 : const casacore::IPosition& pos) const
481 : {
482 0 : ostringstream oss;
483 :
484 0 : oss << "Position = " << pos+1;
485 0 : casacore::String temp(oss);
486 0 : title = temp;
487 0 : }
488 :
489 :
490 : template <class T>
491 0 : void MomentCalcBase<T>::setCoordinateSystem (casacore::CoordinateSystem& cSys,
492 : MomentsBase<T>& iMom)
493 : {
494 0 : cSys = iMom.coordinates() ;
495 0 : }
496 :
497 : template <class T>
498 0 : void MomentCalcBase<T>::setUpCoords (const MomentsBase<T>& iMom,
499 : casacore::Vector<casacore::Double>& pixelIn,
500 : casacore::Vector<casacore::Double>& worldOut,
501 : casacore::Vector<casacore::Double>& sepWorldCoord,
502 : casacore::LogIO& os,
503 : casacore::Double& integratedScaleFactor,
504 : const casacore::CoordinateSystem& cSys,
505 : casacore::Bool doCoordProfile,
506 : casacore::Bool doCoordRandom) const
507 : //
508 : // casacore::Input:
509 : // doCoordProfile - we need the coordinate for each pixel of the profile
510 : // and we precompute it if we can
511 : // doCoordRandom - we need the coordinate for occaisional use
512 : //
513 : // This function does two things. It sets up the pixelIn
514 : // and worldOut vectors needed by getMomentCoord. It also
515 : // precomputes the vector of coordinates for the moment axis
516 : // profile if it is separable
517 : //
518 : {
519 :
520 : // Do we need the scale factor for the integrated moment
521 :
522 0 : casacore::Int axis = iMom.momentAxis_p;
523 0 : casacore::Bool doIntScaleFactor = false;
524 0 : integratedScaleFactor = 1.0;
525 0 : for (casacore::uInt i=0; i<iMom.moments_p.nelements(); i++) {
526 0 : if (iMom.moments_p(i) == MomentsBase<casacore::Float>::INTEGRATED) {
527 0 : doIntScaleFactor = true;
528 0 : break;
529 : }
530 : }
531 : //
532 0 : sepWorldCoord.resize(0);
533 0 : if (!doCoordProfile && !doCoordRandom && !doIntScaleFactor) return;
534 :
535 : // Resize these vectors used for occaisional coordinate transformations
536 :
537 0 : pixelIn.resize(cSys.nPixelAxes());
538 0 : worldOut.resize(cSys.nWorldAxes());
539 0 : if (!doCoordProfile && !doIntScaleFactor) return;
540 :
541 : // Find the coordinate for the moment axis
542 :
543 : casacore::Int coordinate, axisInCoordinate;
544 0 : cSys.findPixelAxis(coordinate, axisInCoordinate, axis);
545 :
546 : // Find out whether this coordinate is separable or not
547 :
548 0 : casacore::Int nPixelAxes = cSys.coordinate(coordinate).nPixelAxes();
549 0 : casacore::Int nWorldAxes = cSys.coordinate(coordinate).nWorldAxes();
550 :
551 : // Precompute the profile coordinates if it is separable and needed
552 : // The Integrated moment scale factor is worked out here as well so the
553 : // logic is a bit contorted
554 :
555 0 : casacore::Bool doneIntScale = false;
556 0 : if (nPixelAxes == 1 && nWorldAxes == 1) {
557 0 : pixelIn = cSys_p.referencePixel();
558 : //
559 0 : casacore::Vector<casacore::Double> frequency(iMom.getShape()(axis));
560 0 : if (doCoordProfile) {
561 0 : for (casacore::uInt i=0; i<frequency.nelements(); i++) {
562 0 : frequency(i) = getMomentCoord(iMom, pixelIn, worldOut, casacore::Double(i));
563 : }
564 : }
565 :
566 : // If the coordinate of the moment axis is Spectral convert to km/s
567 : // Although I could work this out here, it would be decoupled from
568 : // ImageMoments which works the same thing out and sets the units.
569 : // So to ensure coupling, i pass in this switch via the IM object
570 :
571 0 : if (iMom.convertToVelocity_p) {
572 0 : AlwaysAssert(cSys.type(coordinate)==casacore::Coordinate::SPECTRAL, casacore::AipsError); // Should never fail !
573 : //
574 0 : const casacore::SpectralCoordinate& sc = cSys.spectralCoordinate(coordinate);
575 0 : casacore::SpectralCoordinate sc0(sc);
576 :
577 : // Convert
578 :
579 0 : sc0.setVelocity (casacore::String("km/s"), iMom.velocityType_p);
580 0 : if (doCoordProfile) {
581 0 : sc0.frequencyToVelocity (sepWorldCoord, frequency);
582 : }
583 :
584 : // Find increment in world units at reference pixel if needed
585 :
586 0 : if (doIntScaleFactor) {
587 0 : casacore::Quantum<casacore::Double> vel0, vel1;
588 0 : casacore::Double pix0 = sc0.referencePixel()(0) - 0.5;
589 0 : casacore::Double pix1 = sc0.referencePixel()(0) + 0.5;
590 0 : sc0.pixelToVelocity (vel0, pix0);
591 0 : sc0.pixelToVelocity (vel1, pix1);
592 0 : integratedScaleFactor = abs(vel1.getValue() - vel0.getValue());
593 0 : doneIntScale = true;
594 0 : }
595 0 : }
596 0 : } else {
597 : os << casacore::LogIO::NORMAL
598 0 : << "You have asked for a coordinate moment from a non-separable " << endl;
599 0 : os << "axis. This means a coordinate must be computed for each pixel " << endl;
600 0 : os << "of each profile which will cause performance degradation" << casacore::LogIO::POST;
601 : }
602 : //
603 0 : if (doIntScaleFactor && !doneIntScale) {
604 :
605 : // We need the Integrated moment scale factor but the moment
606 : // axis is non-separable
607 :
608 0 : const casacore::Coordinate& c = cSys.coordinate(coordinate);
609 0 : casacore::Double inc = c.increment()(axisInCoordinate);
610 0 : integratedScaleFactor = abs(inc*inc);
611 0 : doneIntScale = true;
612 : }
613 : }
614 :
615 : template <class T>
616 0 : T& MomentCalcBase<T>::stdDeviation(MomentsBase<T>& iMom) const
617 : {
618 0 : return iMom.stdDeviation_p;
619 : }
620 :
621 : // Fill the ouput moments array
622 : template<class T>
623 0 : void MomentCalcBase<T>::setCalcMoments
624 : (const MomentsBase<T>& iMom,
625 : casacore::Vector<T>& calcMoments,
626 : casacore::Vector<casacore::Bool>& calcMomentsMask,
627 : casacore::Vector<casacore::Double>& pixelIn,
628 : casacore::Vector<casacore::Double>& worldOut,
629 : casacore::Bool doCoord,
630 : casacore::Double integratedScaleFactor,
631 : T dMedian,
632 : T vMedian,
633 : casacore::Int nPts,
634 : typename casacore::NumericTraits<T>::PrecisionType s0,
635 : typename casacore::NumericTraits<T>::PrecisionType s1,
636 : typename casacore::NumericTraits<T>::PrecisionType s2,
637 : typename casacore::NumericTraits<T>::PrecisionType s0Sq,
638 : typename casacore::NumericTraits<T>::PrecisionType sumAbsDev,
639 : T dMin,
640 : T dMax,
641 : casacore::Int iMin,
642 : casacore::Int iMax) const
643 : //
644 : // Fill the moments vector
645 : //
646 : // Inputs
647 : // integratedScaleFactor width of a channel in km/s or Hz or whatever
648 : // Outputs:
649 : // calcMoments The moments
650 : //
651 : {
652 : // casacore::Short hand to fish ImageMoments enum values out
653 : // Despite being our friend, we cannot refer to the
654 : // enum values as just, say, "AVERAGE"
655 :
656 : using IM = MomentsBase<casacore::Float>;
657 :
658 : // Normalize and fill moments
659 :
660 0 : calcMomentsMask = true;
661 0 : calcMoments(IM::AVERAGE) = s0 / nPts;
662 0 : calcMoments(IM::INTEGRATED) = s0 * integratedScaleFactor;
663 0 : if (abs(s0) > 0.0) {
664 0 : calcMoments(IM::WEIGHTED_MEAN_COORDINATE) = s1 / s0;
665 : //
666 0 : calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE) =
667 0 : (s2 / s0) - calcMoments(IM::WEIGHTED_MEAN_COORDINATE) *
668 0 : calcMoments(IM::WEIGHTED_MEAN_COORDINATE);
669 0 : calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE) =
670 0 : abs(calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE));
671 0 : if (calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE) > 0.0) {
672 0 : calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE) =
673 0 : sqrt(calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE));
674 : } else {
675 0 : calcMoments(IM::WEIGHTED_DISPERSION_COORDINATE) = 0.0;
676 0 : calcMomentsMask(IM::WEIGHTED_DISPERSION_COORDINATE) = false;
677 : }
678 : } else {
679 0 : calcMomentsMask(IM::WEIGHTED_MEAN_COORDINATE) = false;
680 0 : calcMomentsMask(IM::WEIGHTED_DISPERSION_COORDINATE) = false;
681 : }
682 :
683 : // Standard deviation about mean of I
684 :
685 0 : if (nPts>1 && casacore::Float((s0Sq - s0*s0/nPts)/(nPts-1)) > 0) {
686 0 : calcMoments(IM::STANDARD_DEVIATION) = sqrt((s0Sq - s0*s0/nPts)/(nPts-1));
687 : } else {
688 0 : calcMoments(IM::STANDARD_DEVIATION) = 0;
689 0 : calcMomentsMask(IM::STANDARD_DEVIATION) = false;
690 : }
691 :
692 : // Rms of I
693 :
694 0 : calcMoments(IM::RMS) = sqrt(s0Sq/nPts);
695 :
696 : // Absolute mean deviation
697 :
698 0 : calcMoments(IM::ABS_MEAN_DEVIATION) = sumAbsDev / nPts;
699 :
700 : // Maximum value
701 :
702 0 : calcMoments(IM::MAXIMUM) = dMax;
703 :
704 : // casacore::Coordinate of min/max value
705 :
706 0 : if (doCoord) {
707 0 : calcMoments(IM::MAXIMUM_COORDINATE) = getMomentCoord(
708 : iMom, pixelIn, worldOut, casacore::Double(iMax),
709 0 : iMom.convertToVelocity_p
710 : );
711 0 : calcMoments(IM::MINIMUM_COORDINATE) = getMomentCoord(
712 : iMom, pixelIn, worldOut, casacore::Double(iMin),
713 0 : iMom.convertToVelocity_p
714 : );
715 : } else {
716 0 : calcMoments(IM::MAXIMUM_COORDINATE) = 0.0;
717 0 : calcMoments(IM::MINIMUM_COORDINATE) = 0.0;
718 0 : calcMomentsMask(IM::MAXIMUM_COORDINATE) = false;
719 0 : calcMomentsMask(IM::MINIMUM_COORDINATE) = false;
720 : }
721 :
722 : // Minimum value
723 0 : calcMoments(IM::MINIMUM) = dMin;
724 :
725 : // Medians
726 :
727 0 : calcMoments(IM::MEDIAN) = dMedian;
728 0 : calcMoments(IM::MEDIAN_COORDINATE) = vMedian;
729 0 : }
730 :
731 : }
732 :
|
