Line data Source code
1 : //# MomentWindow.cc:
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 : #include <imageanalysis/ImageAnalysis/MomentWindow.h>
29 :
30 : #include <casacore/casa/aips.h>
31 : #include <casacore/casa/Arrays/Vector.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
34 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
35 : #include <casacore/scimath/Fitting/NonLinearFitLM.h>
36 : #include <casacore/scimath/Functionals/Polynomial.h>
37 : #include <casacore/scimath/Functionals/CompoundFunction.h>
38 : #include <imageanalysis/ImageAnalysis/ImageMoments.h>
39 : #include <casacore/lattices/LatticeMath/LatticeStatsBase.h>
40 : #include <casacore/casa/BasicMath/Math.h>
41 : #include <casacore/casa/Logging/LogIO.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 :
45 : namespace casa {
46 :
47 : // Derived class MomentWindow
48 :
49 : template <class T>
50 0 : MomentWindow<T>::MomentWindow(shared_ptr<casacore::Lattice<T>> pAncilliaryLattice,
51 : MomentsBase<T>& iMom,
52 : casacore::LogIO& os,
53 : const casacore::uInt nLatticeOut)
54 0 : : _ancilliaryLattice(pAncilliaryLattice),
55 0 : iMom_p(iMom),
56 0 : os_p(os)
57 : {
58 : // Set moment selection vector
59 :
60 0 : selectMoments_p = this->selectMoments(iMom_p);
61 :
62 : // Set/check some dimensionality
63 :
64 0 : constructorCheck(calcMoments_p, calcMomentsMask_p, selectMoments_p, nLatticeOut);
65 :
66 : // Fish out moment axis
67 :
68 0 : casacore::Int momAxis = this->momentAxis(iMom_p);
69 :
70 : // Set up slice shape for extraction from masking lattice
71 :
72 0 : if (_ancilliaryLattice != 0) {
73 0 : sliceShape_p.resize(_ancilliaryLattice->ndim());
74 0 : sliceShape_p = 1;
75 0 : sliceShape_p(momAxis) = _ancilliaryLattice->shape()(momAxis);
76 : }
77 :
78 :
79 : // this->yAutoMinMax(yMinAuto_p, yMaxAuto_p, iMom_p);
80 :
81 : // Are we computing the expensive moments ?
82 :
83 0 : this->costlyMoments(iMom_p, doMedianI_p, doMedianV_p, doAbsDev_p);
84 :
85 : // Are we computing coordinate-dependent moments. If
86 : // so precompute coordinate vector is momebt axis separable
87 :
88 0 : this->setCoordinateSystem (cSys_p, iMom_p);
89 0 : this->doCoordCalc(doCoordProfile_p, doCoordRandom_p, iMom_p);
90 0 : this->setUpCoords(iMom_p, pixelIn_p, worldOut_p, sepWorldCoord_p, os_p,
91 0 : integratedScaleFactor_p, cSys_p, doCoordProfile_p,
92 0 : doCoordRandom_p);
93 :
94 : // What is the axis type of the moment axis
95 :
96 0 : momAxisType_p = this->momentAxisName(cSys_p, iMom_p);
97 :
98 : // Are we fitting, automatically or interactively ?
99 :
100 0 : doFit_p = this->doFit(iMom_p);
101 :
102 : // Values to assess if spectrum is all noise or not
103 :
104 0 : peakSNR_p = this->peakSNR(iMom_p);
105 0 : stdDeviation_p = this->stdDeviation(iMom_p);
106 :
107 : // Number of failed Gaussian fits
108 :
109 0 : nFailed_p = 0;
110 0 : }
111 :
112 :
113 : template <class T>
114 0 : MomentWindow<T>::~MomentWindow()
115 0 : {;}
116 :
117 : template <class T>
118 0 : void MomentWindow<T>::process(T&,
119 : casacore::Bool&,
120 : const casacore::Vector<T>&,
121 : const casacore::Vector<casacore::Bool>&,
122 : const casacore::IPosition&)
123 : {
124 0 : throw(casacore::AipsError("MomentWindow<T>::process not implemented"));
125 : }
126 :
127 :
128 : template <class T>
129 0 : void MomentWindow<T>::multiProcess(casacore::Vector<T>& moments,
130 : casacore::Vector<casacore::Bool>& momentsMask,
131 : const casacore::Vector<T>& profileIn,
132 : const casacore::Vector<casacore::Bool>& profileInMask,
133 : const casacore::IPosition& inPos)
134 : //
135 : // Generate windowed moments of this profile.
136 : // The profile comes with its own mask (or a null mask
137 : // which means all good). In addition, we create
138 : // a further mask by applying the clip range to either
139 : // the primary lattice, or the ancilliary lattice (e.g.
140 : // the smoothed lattice)
141 : //
142 : {
143 :
144 : // Fish out the ancilliary image slice if needed. Stupid slice functions
145 : // require me to create the slice empty every time so degenerate
146 : // axes can be chucked out. We set up a pointer to the primary or
147 : // ancilliary vector object that we can use for fast access
148 0 : const T* pProfileSelect = 0;
149 : casacore::Bool deleteIt;
150 0 : if (_ancilliaryLattice) {
151 0 : casacore::Array<T> ancilliarySlice;
152 0 : casacore::IPosition stride(_ancilliaryLattice->ndim(),1);
153 0 : _ancilliaryLattice->getSlice(ancilliarySlice, inPos,
154 0 : sliceShape_p, stride, true);
155 0 : ancilliarySliceRef_p.reference(ancilliarySlice);
156 :
157 0 : pProfileSelect_p = &ancilliarySliceRef_p;
158 0 : pProfileSelect = ancilliarySliceRef_p.getStorage(deleteIt);
159 0 : } else {
160 0 : pProfileSelect_p = &profileIn;
161 0 : pProfileSelect = profileIn.getStorage(deleteIt);
162 : }
163 :
164 :
165 : // Make abcissa and labels
166 :
167 0 : static casacore::Vector<casacore::Int> window(2);
168 : static casacore::Int nPts = 0;
169 :
170 0 : abcissa_p.resize(pProfileSelect_p->size());
171 0 : indgen(abcissa_p);
172 : //this->makeAbcissa (abcissa_p, pProfileSelect_p->nelements());
173 0 : casacore::String xLabel;
174 0 : if (momAxisType_p.empty()) {
175 0 : xLabel = "x (pixels)";
176 : } else {
177 0 : xLabel = momAxisType_p + " (pixels)";
178 : }
179 0 : const casacore::String yLabel("Intensity");
180 0 : casacore::String title;
181 0 : setPosLabel(title, inPos);
182 :
183 :
184 : // Do the window selection
185 :
186 : // Define the window automatically
187 0 : casacore::Vector<T> gaussPars;
188 0 : if (getAutoWindow(nFailed_p, window, abcissa_p, *pProfileSelect_p, profileInMask,
189 0 : peakSNR_p, stdDeviation_p, doFit_p)) {
190 0 : nPts = window(1) - window(0) + 1;
191 : }
192 : else {
193 0 : nPts = 0;
194 : }
195 :
196 :
197 :
198 : // If no points make moments zero and mask
199 :
200 0 : if (nPts==0) {
201 0 : moments = 0.0;
202 0 : momentsMask = false;
203 :
204 0 : if (_ancilliaryLattice) {
205 0 : ancilliarySliceRef_p.freeStorage(pProfileSelect, deleteIt);
206 : } else {
207 0 : profileIn.freeStorage(pProfileSelect, deleteIt);
208 : }
209 0 : return;
210 : }
211 :
212 :
213 : // Resize array for median. Is resized correctly later
214 :
215 0 : selectedData_p.resize(nPts);
216 :
217 :
218 : // Were the profile coordinates precomputed ?
219 :
220 0 : casacore::Bool preComp = (sepWorldCoord_p.nelements() > 0);
221 :
222 : //
223 : // We must fill in the input pixel coordinate if we need
224 : // coordinates, but did not pre compute them
225 : //
226 0 : if (!preComp) {
227 0 : if (doCoordRandom_p || doCoordProfile_p) {
228 0 : for (casacore::uInt i=0; i<inPos.nelements(); i++) {
229 0 : pixelIn_p(i) = casacore::Double(inPos(i));
230 : }
231 : }
232 : }
233 :
234 :
235 : // Set up a pointer for fast access to the profile mask
236 : // if it exists.
237 :
238 : casacore::Bool deleteIt2;
239 0 : const casacore::Bool* pProfileInMask = profileInMask.getStorage(deleteIt2);
240 :
241 :
242 : // Accumulate sums and acquire selected data from primary lattice
243 :
244 0 : typename casacore::NumericTraits<T>::PrecisionType s0 = 0.0;
245 0 : typename casacore::NumericTraits<T>::PrecisionType s0Sq = 0.0;
246 0 : typename casacore::NumericTraits<T>::PrecisionType s1 = 0.0;
247 0 : typename casacore::NumericTraits<T>::PrecisionType s2 = 0.0;
248 0 : casacore::Int iMin = -1;
249 0 : casacore::Int iMax = -1;
250 0 : T dMin = 1.0e30;
251 0 : T dMax = -1.0e30;
252 0 : casacore::Double coord = 0.0;
253 :
254 : casacore::Int i,j;
255 0 : for (i=window(0),j=0; i<=window(1); i++) {
256 0 : if (pProfileInMask[i]) {
257 0 : if (preComp) {
258 0 : coord = sepWorldCoord_p(i);
259 0 : } else if (doCoordProfile_p) {
260 0 : coord = this->getMomentCoord(iMom_p, pixelIn_p,
261 0 : worldOut_p, casacore::Double(i));
262 : }
263 0 : this->accumSums(s0, s0Sq, s1, s2, iMin, iMax,
264 0 : dMin, dMax, i, profileIn(i), coord);
265 0 : selectedData_p(j) = profileIn(i);
266 0 : j++;
267 : }
268 : }
269 0 : nPts = j;
270 :
271 :
272 : // Absolute deviations of I from mean needs an extra pass.
273 :
274 0 : typename casacore::NumericTraits<T>::PrecisionType sumAbsDev = 0.0;
275 0 : if (doAbsDev_p) {
276 0 : T iMean = s0 / nPts;
277 0 : for (i=0; i<nPts; i++) sumAbsDev += abs(selectedData_p(i) - iMean);
278 : }
279 :
280 :
281 :
282 : // Delete memory associated with pointers
283 :
284 0 : if (_ancilliaryLattice) {
285 0 : ancilliarySliceRef_p.freeStorage(pProfileSelect, deleteIt);
286 : } else {
287 0 : profileIn.freeStorage(pProfileSelect, deleteIt);
288 : }
289 0 : profileInMask.freeStorage(pProfileInMask, deleteIt2);
290 :
291 :
292 : // Median of I
293 :
294 0 : T dMedian = 0.0;
295 0 : if (doMedianI_p) {
296 0 : selectedData_p.resize(nPts,true);
297 0 : dMedian = median(selectedData_p);
298 : }
299 :
300 : // Fill all moments array
301 :
302 0 : T vMedian = 0;
303 0 : this->setCalcMoments(iMom_p, calcMoments_p, calcMomentsMask_p, pixelIn_p,
304 0 : worldOut_p, doCoordRandom_p, integratedScaleFactor_p,
305 : dMedian, vMedian, nPts, s0, s1, s2, s0Sq,
306 : sumAbsDev, dMin, dMax, iMin, iMax);
307 :
308 :
309 : // Fill selected moments
310 :
311 0 : for (i=0; i<casacore::Int(selectMoments_p.nelements()); i++) {
312 0 : moments(i) = calcMoments_p(selectMoments_p(i));
313 0 : momentsMask(i) = true;
314 0 : momentsMask(i) = calcMomentsMask_p(selectMoments_p(i));
315 : }
316 0 : }
317 :
318 : template <class T>
319 0 : Bool MomentWindow<T>::getAutoWindow (casacore::uInt& nFailed,
320 : casacore::Vector<casacore::Int>& window,
321 : const casacore::Vector<T>& x,
322 : const casacore::Vector<T>& y,
323 : const casacore::Vector<casacore::Bool>& mask,
324 : const T peakSNR,
325 : const T stdDeviation,
326 : const casacore::Bool doFit) const
327 : //
328 : // Automatically fit a Gaussian and return the +/- 3-sigma window or
329 : // invoke Bosma's method to set a window. If a plotting device is
330 : // active, we also plot the spectra and fits
331 : //
332 : // Inputs:
333 : // x,y Spectrum
334 : // mask Mask associated with spectrum. true is good.
335 : // plotter Plot spectrum and optionally the window
336 : // x,yLabel x label for plots
337 : // title
338 : // casacore::Input/output
339 : // nFailed Cumulative number of failed fits
340 : // Output:
341 : // window The window (pixels). If both 0, then discard this spectrum
342 : // and mask moments
343 : //
344 : {
345 0 : if (doFit) {
346 0 : casacore::Vector<T> gaussPars(4);
347 0 : if (!this->getAutoGaussianFit (nFailed, gaussPars, x, y, mask, peakSNR, stdDeviation)) {
348 0 : window = 0;
349 0 : return false;
350 : } else {
351 :
352 : // Set 3-sigma limits. This assumes that there are some unmasked
353 : // points in the window !
354 :
355 0 : if (!setNSigmaWindow (window, gaussPars(1), gaussPars(2),
356 0 : y.nelements(), 3)) {
357 0 : window = 0;
358 0 : return false;
359 : }
360 : }
361 0 : } else {
362 : // Invoke Albert's method (see AJ, 86, 1791)
363 :
364 0 : if (! _getBosmaWindow (window, y, mask, peakSNR, stdDeviation)) {
365 0 : window = 0;
366 0 : return false;
367 : }
368 : }
369 0 : return true;
370 : }
371 :
372 0 : template <class T> casacore::Bool MomentWindow<T>::_getBosmaWindow(
373 : casacore::Vector<casacore::Int>& window, const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask,
374 : const T peakSNR, const T stdDeviation
375 : ) const {
376 : // Automatically work out the spectral window
377 : // with Albert Bosma's algorithm.
378 : // Inputs:
379 : // x,y Spectrum
380 : // Output:
381 : // window The window
382 : // casacore::Bool false if we reject this spectrum. This may
383 : // be because it is all noise, or all masked
384 :
385 : // See if this spectrum is all noise first. If so, forget it.
386 : // Return straight away if all maske
387 : T dMean;
388 0 : casacore::uInt iNoise = this->allNoise(dMean, y, mask, peakSNR, stdDeviation);
389 0 : if (iNoise == 2) {
390 : // all masked
391 0 : return false;
392 : }
393 :
394 0 : if (iNoise==1) {
395 : // all noise
396 0 : window = 0;
397 0 : return false;
398 : }
399 : // Find peak
400 :
401 0 : casacore::ClassicalStatistics<AccumType, DataIterator, MaskIterator> statsCalculator;
402 0 : statsCalculator.addData(y.begin(), mask.begin(), y.size());
403 0 : StatsData<AccumType> stats = statsCalculator.getStatistics();
404 0 : const casacore::Int nPts = y.size();
405 0 : auto maxPos = stats.maxpos.second;
406 0 : casacore::Int iMin = max(0, casacore::Int(maxPos)-2);
407 0 : casacore::Int iMax = min(nPts-1, casacore::Int(maxPos)+2);
408 0 : T tol = stdDeviation / (nPts - (iMax-iMin-1));
409 : // Iterate to convergence
410 0 : auto first = true;
411 0 : auto converged = false;
412 0 : auto more = true;
413 0 : T yMean = 0;
414 0 : T oldYMean = 0;
415 0 : while (more) {
416 : // Find mean outside of peak region
417 0 : AccumType sum = 0;
418 : casacore::Int i,j;
419 0 : for (i=0,j=0; i<nPts; ++i) {
420 0 : if (mask[i] && (i < iMin || i > iMax)) {
421 0 : sum += y[i];
422 0 : ++j;
423 : }
424 : }
425 0 : if (j>0) {
426 0 : yMean = sum / j;
427 : }
428 : // Interpret result
429 0 : if (!first && j>0 && abs(yMean-oldYMean) < tol) {
430 0 : converged = true;
431 0 : more = false;
432 : }
433 0 : else if (iMin==0 && iMax==nPts-1) {
434 0 : more = false;
435 : }
436 : else {
437 : // Widen window and redetermine tolerance
438 0 : oldYMean = yMean;
439 0 : iMin = max(0,iMin - 2);
440 0 : iMax = min(nPts-1,iMax+2);
441 0 : tol = stdDeviation / (nPts - (iMax-iMin-1));
442 : }
443 0 : first = false;
444 : }
445 :
446 : // Return window
447 :
448 0 : if (converged) {
449 0 : window[0] = iMin;
450 0 : window[1] = iMax;
451 0 : return true;
452 : }
453 : else {
454 0 : window = 0;
455 0 : return false;
456 : }
457 0 : }
458 :
459 : template <class T>
460 0 : Bool MomentWindow<T>::setNSigmaWindow (casacore::Vector<casacore::Int>& window,
461 : const T pos,
462 : const T width,
463 : const casacore::Int nPts,
464 : const casacore::Int N) const
465 : //
466 : // Take the fitted Gaussian position and width and
467 : // set an N-sigma window. If the window is too small
468 : // return a Fail condition.
469 : //
470 : // Inputs:
471 : // pos,width The position and width in pixels
472 : // nPts The number of points in the spectrum that was fit
473 : // N The N-sigma
474 : // Outputs:
475 : // window The window in pixels
476 : // casacore::Bool false if window too small to be sensible
477 : //
478 : {
479 0 : window(0) = casacore::Int((pos-N*width)+0.5);
480 0 : window(0) = min(nPts-1,max(0,window(0)));
481 0 : window(1) = casacore::Int((pos+N*width)+0.5);
482 0 : window(1) = min(nPts-1,max(0,window(1)));
483 : // FIXME this was
484 : // if ( abs(window(1)-window(0)) < 3) return false;
485 : // return true;
486 : // but because window(1) - window(0) could be negative and true could be
487 : // returned, an allocation error was occuring because in another function a
488 : // vector was being resized to (window(1) - window(0)). It is possible that
489 : // in that case the absolute value should be calculated but I don't have time
490 : // at the moment to trace through the code and make sure that is really the
491 : // correct thing to do. Thus, making this function return false if window(1) - window(0)
492 : // seems the more conservative approach, so I'm doing that for now.
493 0 : return window(1)-window(0) >= 3;
494 :
495 : }
496 :
497 : }
|
