Line data Source code
1 : //# MomentClip.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: 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.tcc 19940 2007-02-27 05:35:22Z Malte.Marquarding $
27 : //
28 :
29 : #include <imageanalysis/ImageAnalysis/MomentClip.h>
30 :
31 : namespace casa {
32 :
33 0 : template <class T> MomentClip<T>::MomentClip(
34 : shared_ptr<casacore::Lattice<T>> pAncilliaryLattice, MomentsBase<T>& iMom,
35 : casacore::LogIO& os, const casacore::uInt nLatticeOut
36 0 : ) : _ancilliaryLattice(pAncilliaryLattice), iMom_p(iMom), os_p(os) {
37 0 : selectMoments_p = this->selectMoments(iMom_p);
38 0 : constructorCheck(calcMoments_p, calcMomentsMask_p, selectMoments_p, nLatticeOut);
39 0 : casacore::Int momAxis = this->momentAxis(iMom_p);
40 0 : if (_ancilliaryLattice) {
41 0 : sliceShape_p.resize(_ancilliaryLattice->ndim());
42 0 : sliceShape_p = 1;
43 0 : sliceShape_p(momAxis) = _ancilliaryLattice->shape()(momAxis);
44 : }
45 : //this->yAutoMinMax(yMinAuto_p, yMaxAuto_p, iMom_p);
46 :
47 0 : this->selectRange(range_p, doInclude_p, doExclude_p, iMom_p);
48 :
49 0 : this->costlyMoments(iMom_p, doMedianI_p, doMedianV_p, doAbsDev_p);
50 :
51 : // Are we computing coordinate-dependent moments. If so
52 : // precompute coordinate vector if moment axis separable
53 :
54 0 : this->setCoordinateSystem (cSys_p, iMom_p);
55 0 : this->doCoordCalc(doCoordProfile_p, doCoordRandom_p, iMom_p);
56 0 : this->setUpCoords(
57 0 : iMom_p, pixelIn_p, worldOut_p, sepWorldCoord_p, os_p,
58 0 : integratedScaleFactor_p, cSys_p, doCoordProfile_p,
59 0 : doCoordRandom_p
60 : );
61 :
62 : // What is the axis type of the moment axis
63 :
64 0 : momAxisType_p = this->momentAxisName(cSys_p, iMom_p);
65 :
66 : // Number of failed Gaussian fits
67 0 : nFailed_p = 0;
68 0 : }
69 :
70 0 : template <class T> MomentClip<T>::~MomentClip() {}
71 :
72 0 : template <class T> void MomentClip<T>::process(
73 : T&, casacore::Bool&, const casacore::Vector<T>&, const casacore::Vector<casacore::Bool>&,
74 : const casacore::IPosition&
75 : ) {
76 0 : ThrowCc("MomentClip<T>::process(Vector<T>&, IPosition&): not implemented");
77 : }
78 :
79 0 : template <class T> void MomentClip<T>::multiProcess(
80 : casacore::Vector<T>& moments, casacore::Vector<casacore::Bool>& momentsMask,
81 : const casacore::Vector<T>& profileIn, const casacore::Vector<casacore::Bool>& profileInMask,
82 : const casacore::IPosition& inPos
83 : ) {
84 : // The profile comes with its own mask (or a null mask
85 : // which means all good). In addition, we create
86 : // a further mask by applying the clip range to either
87 : // the primary lattice, or the ancilliary lattice (e.g.
88 : // the smoothed lattice)
89 :
90 : // Fish out the ancilliary image slice if needed. Stupid slice functions
91 : // require me to create the slice empty every time so degenerate
92 : // axes can be chucked out. We set up a pointer to the primary or
93 : // ancilliary vector object that we can use for fast access
94 :
95 0 : const T* pProfileSelect = nullptr;
96 0 : auto deleteIt = false;
97 0 : if (_ancilliaryLattice && (doInclude_p || doExclude_p)) {
98 0 : casacore::Array<T> ancilliarySlice;
99 0 : casacore::IPosition stride(_ancilliaryLattice->ndim(),1);
100 0 : _ancilliaryLattice->getSlice(
101 0 : ancilliarySlice, inPos, sliceShape_p, stride, true
102 : );
103 0 : ancilliarySliceRef_p.reference(ancilliarySlice);
104 0 : pProfileSelect_p = &ancilliarySliceRef_p;
105 0 : pProfileSelect = ancilliarySliceRef_p.getStorage(deleteIt);
106 : }
107 : else {
108 0 : pProfileSelect_p = &profileIn;
109 0 : pProfileSelect = profileIn.getStorage(deleteIt);
110 : }
111 : // Resize array for median. Is resized correctly later
112 0 : auto nPts = profileIn.size();
113 0 : selectedData_p.resize(nPts);
114 0 : selectedDataIndex_p.resize(nPts);
115 :
116 : // Were the profile coordinates precomputed ?
117 :
118 0 : auto preComp = sepWorldCoord_p.nelements() > 0;
119 :
120 : // We must fill in the input pixel coordinate if we need
121 : // coordinates, but did not pre compute them
122 0 : if (! preComp && (doCoordRandom_p || doCoordProfile_p)) {
123 0 : for (casacore::uInt i=0; i<inPos.size(); ++i) {
124 0 : pixelIn_p[i] = casacore::Double(inPos[i]);
125 : }
126 : }
127 : // Compute moments. The actual moment computation always done with
128 : // the original data, regardless of whether the pixel selection is
129 : // done with the primary or ancilliary data.
130 :
131 0 : typename casacore::NumericTraits<T>::PrecisionType s0 = 0.0;
132 0 : typename casacore::NumericTraits<T>::PrecisionType s0Sq = 0.0;
133 0 : typename casacore::NumericTraits<T>::PrecisionType s1 = 0.0;
134 0 : typename casacore::NumericTraits<T>::PrecisionType s2 = 0.0;
135 0 : casacore::Int iMin = -1;
136 0 : casacore::Int iMax = -1;
137 0 : T dMin = 1.0e30;
138 0 : T dMax = -1.0e30;
139 0 : casacore::Double coord = 0.0;
140 : casacore::uInt i, j;
141 0 : if (profileInMask.empty()) {
142 : // No mask included.
143 0 : if (doInclude_p) {
144 0 : for (i=0,j=0; i<nPts; ++i) {
145 0 : if (
146 0 : pProfileSelect[i] >= range_p[0]
147 0 : && pProfileSelect[i] <= range_p[1]
148 : ) {
149 0 : if (preComp) {
150 0 : coord = sepWorldCoord_p(i);
151 : }
152 0 : else if (doCoordProfile_p) {
153 0 : coord = this->getMomentCoord(
154 0 : iMom_p, pixelIn_p,
155 0 : worldOut_p, casacore::Double(i)
156 : );
157 : }
158 0 : this->accumSums(
159 : s0, s0Sq, s1, s2, iMin, iMax,
160 0 : dMin, dMax, i, profileIn(i), coord
161 : );
162 0 : selectedData_p[j] = profileIn[i];
163 0 : selectedDataIndex_p[j] = i;
164 0 : ++j;
165 : }
166 : }
167 0 : nPts = j;
168 : }
169 0 : else if (doExclude_p) {
170 0 : for (i=0,j=0; i<nPts; ++i) {
171 0 : if (
172 0 : pProfileSelect[i] <= range_p[0]
173 0 : || pProfileSelect[i] >= range_p[1]
174 : ) {
175 0 : if (preComp) {
176 0 : coord = sepWorldCoord_p[i];
177 : }
178 0 : else if (doCoordProfile_p) {
179 0 : coord = this->getMomentCoord(
180 0 : iMom_p, pixelIn_p,
181 0 : worldOut_p, casacore::Double(i)
182 : );
183 : }
184 0 : this->accumSums(
185 : s0, s0Sq, s1, s2, iMin, iMax,
186 0 : dMin, dMax, i, profileIn[i], coord
187 : );
188 0 : selectedData_p[j] = profileIn[i];
189 0 : selectedDataIndex_p[j] = i;
190 0 : ++j;
191 : }
192 : }
193 0 : nPts = j;
194 : }
195 : else {
196 0 : for (i=0; i<nPts; ++i) {
197 0 : if (preComp) {
198 0 : coord = sepWorldCoord_p[i];
199 : }
200 0 : else if (doCoordProfile_p) {
201 0 : coord = this->getMomentCoord(
202 0 : iMom_p, pixelIn_p,
203 0 : worldOut_p, casacore::Double(i)
204 : );
205 : }
206 0 : this->accumSums(
207 : s0, s0Sq, s1, s2, iMin, iMax,
208 0 : dMin, dMax, i, profileIn[i], coord
209 : );
210 0 : selectedData_p[i] = profileIn[i];
211 0 : selectedDataIndex_p[i] = i;
212 : }
213 : }
214 : }
215 : else {
216 : // Set up a pointer for faster access to the profile mask
217 0 : auto deleteIt2 = false;
218 0 : const auto* pProfileInMask = profileInMask.getStorage(deleteIt2);
219 :
220 0 : if (doInclude_p) {
221 0 : for (i=0, j=0; i<nPts; ++i) {
222 0 : if (
223 : pProfileInMask[i]
224 0 : && pProfileSelect[i] >= range_p(0)
225 0 : && pProfileSelect[i] <= range_p(1)
226 : ) {
227 0 : if (preComp) {
228 0 : coord = sepWorldCoord_p[i];
229 : }
230 0 : else if (doCoordProfile_p) {
231 0 : coord = this->getMomentCoord(
232 0 : iMom_p, pixelIn_p, worldOut_p, casacore::Double(i)
233 : );
234 : }
235 0 : this->accumSums(
236 : s0, s0Sq, s1, s2, iMin, iMax,
237 0 : dMin, dMax, i, profileIn(i), coord
238 : );
239 0 : selectedData_p[j] = profileIn[i];
240 0 : selectedDataIndex_p[j] = i;
241 0 : ++j;
242 : }
243 : }
244 : }
245 0 : else if (doExclude_p) {
246 0 : for (i=0, j=0; i<nPts; ++i) {
247 0 : if (
248 : pProfileInMask[i]
249 0 : && (
250 0 : pProfileSelect[i] <= range_p[0]
251 0 : || pProfileSelect[i] >= range_p[1]
252 : )
253 : ) {
254 0 : if (preComp) {
255 0 : coord = sepWorldCoord_p(i);
256 : }
257 0 : else if (doCoordProfile_p) {
258 0 : coord = this->getMomentCoord(
259 0 : iMom_p, pixelIn_p,
260 0 : worldOut_p, casacore::Double(i)
261 : );
262 : }
263 0 : this->accumSums(
264 : s0, s0Sq, s1, s2, iMin, iMax,
265 0 : dMin, dMax, i, profileIn[i], coord
266 : );
267 0 : selectedData_p[j] = profileIn[i];
268 0 : selectedDataIndex_p[j] = i;
269 0 : ++j;
270 : }
271 : }
272 : } else {
273 0 : for (i=0,j=0; i<nPts; ++i) {
274 0 : if (pProfileInMask[i]) {
275 0 : if (preComp) {
276 0 : coord = sepWorldCoord_p[i];
277 : }
278 0 : else if (doCoordProfile_p) {
279 0 : coord = this->getMomentCoord(
280 0 : iMom_p, pixelIn_p,
281 0 : worldOut_p, casacore::Double(i)
282 : );
283 : }
284 0 : this->accumSums(
285 : s0, s0Sq, s1, s2, iMin, iMax,
286 0 : dMin, dMax, i, profileIn[i], coord
287 : );
288 0 : selectedData_p[j] = profileIn[i];
289 0 : selectedDataIndex_p[j] = i;
290 0 : ++j;
291 : }
292 : }
293 : }
294 0 : nPts = j;
295 0 : profileInMask.freeStorage(pProfileInMask, deleteIt2);
296 : }
297 : // Delete pointer memory
298 0 : if (_ancilliaryLattice && (doInclude_p || doExclude_p)) {
299 0 : ancilliarySliceRef_p.freeStorage(pProfileSelect, deleteIt);
300 : }
301 : else {
302 0 : profileIn.freeStorage(pProfileSelect, deleteIt);
303 : }
304 : // If no points make moments zero and mask
305 0 : if (nPts == 0) {
306 0 : moments = 0.0;
307 0 : momentsMask = false;
308 0 : return;
309 : }
310 : // Absolute deviations of I from mean needs an extra pass.
311 0 : typename casacore::NumericTraits<T>::PrecisionType sumAbsDev = 0;
312 0 : if (doAbsDev_p) {
313 0 : T iMean = s0 / nPts;
314 0 : for (i=0; i<nPts; ++i) {
315 0 : sumAbsDev += abs(selectedData_p(i) - iMean);
316 : }
317 : }
318 : // Median of I
319 0 : T dMedian = 0.0;
320 0 : if (doMedianI_p) {
321 0 : selectedData_p.resize(nPts,true);
322 0 : dMedian = median(selectedData_p);
323 : }
324 : // Median coordinate. ImageMoments will only be allowing this if
325 : // we are not offering the ancilliary lattice, and with an include
326 : // or exclude range. Pretty dodgy
327 0 : T vMedian(0.0);
328 0 : if (doMedianV_p) {
329 0 : if (doInclude_p || doExclude_p) {
330 : // Treat spectrum as a probability distribution for velocity
331 : // and generate cumulative probability (it's already sorted
332 : // of course).
333 0 : selectedData_p.resize(nPts,true);
334 0 : selectedData_p(0) = abs(selectedData_p[0]);
335 0 : auto dataMax = selectedData_p[0];
336 0 : for (i=1; i<nPts; ++i) {
337 0 : selectedData_p[i] += abs(selectedData_p[i-1]);
338 0 : dataMax = max(dataMax, selectedData_p[i]);
339 : }
340 : // Find 1/2 way value (well, the first one that occurs)
341 :
342 0 : auto halfMax = dataMax/2.0;
343 0 : casacore::Int iVal = 0;
344 0 : for (i=0; i<nPts; ++i) {
345 0 : if (selectedData_p[i] >= halfMax) {
346 0 : iVal = i;
347 0 : break;
348 : }
349 : }
350 : // Linearly interpolate to velocity index
351 : casacore::Double interpPixel;
352 0 : if (iVal > 0) {
353 0 : casacore::Double m = (selectedData_p[iVal] - selectedData_p[iVal-1])
354 0 : / (selectedDataIndex_p[iVal] - selectedDataIndex_p[iVal-1]);
355 0 : casacore::Double b = selectedData_p[iVal] - m*selectedDataIndex_p[iVal];
356 0 : interpPixel = (selectedData_p[iVal] - b)/m;
357 : }
358 : else {
359 0 : interpPixel = selectedDataIndex_p[iVal];
360 : }
361 : // Find world coordinate of that pixel on the moment axis
362 0 : vMedian = this->getMomentCoord(
363 0 : iMom_p, pixelIn_p, worldOut_p, interpPixel,
364 0 : iMom_p.shouldConvertToVelocity()
365 : );
366 : }
367 : }
368 : // Fill all moments array
369 0 : this->setCalcMoments(
370 0 : iMom_p, calcMoments_p, calcMomentsMask_p, pixelIn_p, worldOut_p,
371 0 : doCoordRandom_p, integratedScaleFactor_p,
372 : dMedian, vMedian, nPts, s0, s1, s2, s0Sq,
373 : sumAbsDev, dMin, dMax, iMin, iMax
374 : );
375 : // Fill vector of selected moments
376 0 : for (i=0; i<selectMoments_p.size(); ++i) {
377 0 : moments[i] = calcMoments_p(selectMoments_p[i]);
378 0 : momentsMask[i] = calcMomentsMask_p(selectMoments_p[i]);
379 : }
380 : }
381 :
382 : }
383 :
|
