Line data Source code
1 : //# MomentFit.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 : #include <imageanalysis/ImageAnalysis/MomentFit.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 MomentFit
48 :
49 : template <class T>
50 0 : MomentFit<T>::MomentFit(
51 : MomentsBase<T>& iMom,
52 : casacore::LogIO& os,
53 : const casacore::uInt nLatticeOut)
54 : : iMom_p(iMom),
55 0 : os_p(os)
56 : {
57 : // Set moment selection vector
58 :
59 0 : selectMoments_p = this->selectMoments(iMom_p);
60 :
61 : // Set/check some dimensionality
62 :
63 0 : constructorCheck(calcMoments_p, calcMomentsMask_p, selectMoments_p, nLatticeOut);
64 :
65 : //this->yAutoMinMax(yMinAuto_p, yMaxAuto_p, iMom_p);
66 :
67 : // Are we computing the expensive moments ?
68 :
69 0 : this->costlyMoments(iMom_p, doMedianI_p, doMedianV_p, doAbsDev_p);
70 :
71 : // Are we computing coordinate-dependent moments. If so
72 : // precompute coordinate vector if moment axis is separable
73 :
74 0 : this->setCoordinateSystem (cSys_p, iMom_p);
75 0 : this->doCoordCalc(doCoordProfile_p, doCoordRandom_p, iMom_p);
76 0 : this->setUpCoords(iMom_p, pixelIn_p, worldOut_p, sepWorldCoord_p, os_p,
77 0 : integratedScaleFactor_p, cSys_p, doCoordProfile_p,
78 0 : doCoordRandom_p);
79 :
80 : // What is the axis type of the moment axis
81 :
82 0 : momAxisType_p = this->momentAxisName(cSys_p, iMom_p);
83 :
84 : // Are we fitting, automatically or interactively ?
85 :
86 0 : doFit_p = this->doFit(iMom_p);
87 :
88 : // Values to assess if spectrum is all noise or not
89 :
90 0 : peakSNR_p = this->peakSNR(iMom_p);
91 0 : stdDeviation_p = this->stdDeviation(iMom_p);
92 :
93 : // Number of failed Gaussian fits
94 :
95 0 : nFailed_p = 0;
96 :
97 0 : }
98 :
99 :
100 : template <class T>
101 0 : MomentFit<T>::~MomentFit()
102 0 : {;}
103 :
104 :
105 : template <class T>
106 0 : void MomentFit<T>::process(T&,
107 : casacore::Bool&,
108 : const casacore::Vector<T>&,
109 : const casacore::Vector<casacore::Bool>&,
110 : const casacore::IPosition&)
111 : {
112 0 : throw(casacore::AipsError("MomentFit<T>::process not implemented"));
113 : }
114 :
115 0 : template <class T> void MomentFit<T>::multiProcess(
116 : casacore::Vector<T>& moments,
117 : casacore::Vector<casacore::Bool>& momentsMask,
118 : const casacore::Vector<T>& profileIn,
119 : const casacore::Vector<casacore::Bool>& profileInMask,
120 : const casacore::IPosition& inPos
121 : ) {
122 : // Generate moments from a Gaussian fit of this profile
123 :
124 0 : auto nPts = profileIn.size();
125 0 : casacore::Vector<T> gaussPars(4);
126 0 : abcissa_p.resize(nPts);
127 0 : indgen(abcissa_p);
128 0 : if (
129 0 : ! this->getAutoGaussianFit (
130 0 : nFailed_p, gaussPars, abcissa_p, profileIn, profileInMask,
131 : peakSNR_p, stdDeviation_p
132 : )
133 : ) {
134 0 : moments = 0;
135 0 : momentsMask = false;
136 0 : return;
137 : }
138 : // Were the profile coordinates precomputed ?
139 0 : auto preComp = sepWorldCoord_p.size() > 0;
140 : //
141 : // We must fill in the input pixel coordinate if we need
142 : // coordinates, but did not pre compute them
143 : //
144 0 : if (! preComp && (doCoordRandom_p || doCoordProfile_p)) {
145 0 : for (casacore::uInt i=0; i<inPos.size(); ++i) {
146 0 : pixelIn_p[i] = casacore::Double(inPos[i]);
147 : }
148 : }
149 : // Set Gaussian functional values. We reuse the same functional that
150 : // was used in the interactive fitting display process.
151 0 : gauss_p.setHeight(gaussPars(0));
152 0 : gauss_p.setCenter(gaussPars(1));
153 0 : gauss_p.setWidth(gaussPars(2));
154 :
155 : // Compute moments from the fitted Gaussian
156 :
157 0 : typename casacore::NumericTraits<T>::PrecisionType s0 = 0.0;
158 0 : typename casacore::NumericTraits<T>::PrecisionType s0Sq = 0.0;
159 0 : typename casacore::NumericTraits<T>::PrecisionType s1 = 0.0;
160 0 : typename casacore::NumericTraits<T>::PrecisionType s2 = 0.0;
161 0 : casacore::Int iMin = -1;
162 0 : casacore::Int iMax = -1;
163 0 : T dMin = 1.0e30;
164 0 : T dMax = -1.0e30;
165 0 : casacore::Double coord = 0.0;
166 : T xx;
167 0 : casacore::Vector<T> gData(nPts);
168 :
169 : casacore::uInt i,j;
170 0 : for (i=0, j=0; i<nPts; ++i) {
171 0 : if (profileInMask(i)) {
172 0 : xx = i;
173 0 : gData[j] = gauss_p(xx) + gaussPars[3];
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, gData(j), coord
187 : );
188 0 : ++j;
189 : }
190 : }
191 : // If no unmasked points go home. This shouldn't happen
192 : // as we can't have done a fit under these conditions.
193 0 : nPts = j;
194 0 : if (nPts == 0) {
195 0 : moments = 0;
196 0 : momentsMask = false;
197 0 : return;
198 : }
199 : // Absolute deviations of I from mean needs an extra pass.
200 0 : typename casacore::NumericTraits<T>::PrecisionType sumAbsDev = 0.0;
201 0 : if (doAbsDev_p) {
202 0 : T iMean = s0 / nPts;
203 0 : for (i=0; i<nPts; ++i) {
204 0 : sumAbsDev += abs(gData[i] - iMean);
205 : }
206 : }
207 : // Median of I
208 0 : T dMedian = 0.0;
209 0 : if (doMedianI_p) {
210 0 : gData.resize(nPts, true);
211 0 : dMedian = median(gData);
212 : }
213 0 : T vMedian = 0.0;
214 :
215 : // Fill all moments array
216 :
217 0 : this->setCalcMoments(
218 0 : iMom_p, calcMoments_p, calcMomentsMask_p, pixelIn_p,
219 0 : worldOut_p, doCoordRandom_p, integratedScaleFactor_p,
220 : dMedian, vMedian, nPts, s0, s1, s2, s0Sq,
221 : sumAbsDev, dMin, dMax, iMin, iMax
222 : );
223 :
224 : // Fill vector of selected moments
225 :
226 0 : for (i=0; i<selectMoments_p.size(); ++i) {
227 0 : moments[i] = calcMoments_p(selectMoments_p[i]);
228 0 : momentsMask[i] = true;
229 0 : momentsMask[i] = calcMomentsMask_p(selectMoments_p[i]);
230 : }
231 : }
232 :
233 : }
|
