source: trunk/external-alma/components/SpectralComponents/Spectral2Estimate.tcc@ 3087

Last change on this file since 3087 was 2980, checked in by Malte Marquarding, 10 years ago

Add a copy of casacore/components/SpectralComponents to external-alma directory to prepare for its removal from casacore-trunk. DOn;t activate in SConscript yet.
File size: 8.7 KB
Line 
1//# Spectral2Estimate.cc: Member templates for SpectralEstimate
2//# Copyright (C) 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: Spectral2Estimate.tcc 21465 2014-06-19 05:56:56Z gervandiepen $
27
28//# Includes
29#include <components/SpectralComponents/SpectralEstimate.h>
30
31#include <casa/BasicMath/Math.h>
32#include <casa/BasicSL/Constants.h>
33#include <casa/Utilities/Assert.h>
34#include <components/SpectralComponents/CompiledSpectralElement.h>
35#include <components/SpectralComponents/GaussianSpectralElement.h>
36#include <components/SpectralComponents/PolynomialSpectralElement.h>
37
38namespace casa { //#Begin namespace casa
39
40//# Member templates
41template <class MT>
42const SpectralList &SpectralEstimate::estimate(const Vector<MT> &prof,
43 Vector<MT> *der) {
44 if (prof.nelements() != lprof_p) {
45 delete [] deriv_p; deriv_p = 0; lprof_p = 0;
46 lprof_p = prof.nelements();
47 deriv_p = new Double[lprof_p];
48 };
49 // Check if signal in window
50 if (!window(prof)) return slist_p;
51 // Limit window
52 windowEnd_p = min(windowEnd_p+q_p , Int(lprof_p));
53 windowLow_p = max(windowLow_p-q_p , 0 );
54 // Get the second derivatives
55 findc2(prof);
56 // Next for debugging
57 if (der) {
58 for (uInt i=0; i<lprof_p; i++) (*der)[i] = deriv_p[i];
59 };
60 // Find the estimates (sorted)
61 findga(prof);
62 // cout << slist_p << endl;
63 return slist_p;
64}
65
66template <class MT>
67const SpectralList& SpectralEstimate::estimate(const Vector<MT>& x,
68 const Vector<MT>& y)
69{
70 if (x.nelements() != y.nelements()) {
71 throw(AipsError("Abcissa and ordinate vectors must be the same length"));
72 }
73 if (x.nelements()==1) {
74 throw(AipsError("Not enough elements in vectors"));
75 }
76 // Get pixel-based estimate (into slist_p)
77 estimate(y);
78 // Convert
79 for (uInt i=0; i<slist_p.nelements(); i++) {
80 if (slist_p[i]->getType() != SpectralElement::GAUSSIAN) {
81 throw AipsError("Non-gaussian spectral types cannot be estimated");
82 }
83 const GaussianSpectralElement elIn = *dynamic_cast<const GaussianSpectralElement *>(slist_p[i]);
84 GaussianSpectralElement elOut = convertElement (x, elIn);
85 slist_p.set(elOut, i);
86 }
87 return slist_p;
88}
89
90
91template <class MT>
92uInt SpectralEstimate::window(const Vector<MT> &prof) {
93 windowLow_p =0;
94 windowEnd_p = 0;
95 if (!useWindow_p || rms_p <= 0.0 || lprof_p == 0) {
96 if (regionEnd_p) {
97 windowLow_p = min(max(0,regionLow_p),Int(lprof_p));
98 windowEnd_p = min(regionEnd_p, Int(lprof_p));
99 } else windowEnd_p = lprof_p;
100 return windowEnd_p-windowLow_p;
101 };
102 // Total flux in profile and max position
103 Double flux(0.0);
104 Double pmax(prof(0));
105 uInt imax(0);
106 for (Int i=windowLow_p; i<windowEnd_p; i++) {
107 if (prof(i)>pmax) {
108 pmax = prof(i);
109 imax = i;
110 };
111 flux += prof(i);
112 };
113 // No data
114 if (pmax < cutoff_p) return 0;
115 // Window boundaries; new/old base and centre; width
116 Int width(-1);
117 Int nw(0);
118 Double bnew(flux), bold;
119 Double cnew(imax), cold;
120 do {
121 width++;
122 cold = cnew;
123 bold = bnew;
124 windowLow_p = max(0, Int(cold-width+0.5));
125 windowEnd_p = min(Int(lprof_p), Int(cold+width+1.5));
126 // flux and first moment in window
127 Double s(0);
128 Double c(0);
129 for (Int i=windowLow_p; i<windowEnd_p; i++) {
130 s += prof(i);
131 c += i*prof(i);
132 };
133 bnew = flux-s;
134 nw = lprof_p-windowEnd_p+windowLow_p;
135 if (s != 0.0) {
136 cnew = c/s;
137 if (cnew < 0 || cnew >= lprof_p) cnew = cold;
138 };
139 } while (abs(bnew-bold) > rms_p && nw);
140 return windowEnd_p-windowLow_p;
141}
142
143template <class MT>
144void SpectralEstimate::findc2(const Vector<MT> &prof) {
145 for (Int i=windowLow_p; i<windowEnd_p; i++) {
146 // Moments
147 Double m0(0.0);
148 Double m2(0.0);
149 for (Int j = -q_p; j <= q_p; j++) {
150 Int k = i+j;
151 if (k >= 0 && k<Int(lprof_p)) {
152 // add to moments
153 m0 += prof(k);
154 m2 += prof(k)*j*j;
155 };
156 };
157 // get the derivative
158 deriv_p[i] = a_p*(m2-b_p*m0);
159 };
160}
161
162template <class MT>
163void SpectralEstimate::findga(const Vector<MT> &prof) {
164 Int i(windowLow_p-1);
165 // Window on Gaussian
166 Int iclo(windowLow_p);
167 Int ichi(windowLow_p);
168 // Peak counter
169 Int nmax = 0;
170 GaussianSpectralElement tspel;
171 while (++i < windowEnd_p) {
172 if (deriv_p[i] > 0.0) {
173 // At edge?
174 if (i > windowLow_p && i < windowEnd_p-1) {
175 // Peak in 2nd derivative
176 if (deriv_p[i-1] < deriv_p[i] && deriv_p[i+1] < deriv_p[i]) nmax++;
177 // At start
178 } else if (i == windowLow_p && deriv_p[i+1] < deriv_p[i]) nmax++;
179 // At end of window
180 else if (i == windowEnd_p-1 && deriv_p[i-1] < deriv_p[i]) nmax++;
181 };
182 switch (nmax) {
183 // Search for next peak
184 case 1:
185 break;
186 // Found a Gaussian
187 case 2: {
188 // Some moments
189 Double m0m(0);
190 Double m0(0);
191 Double m1(0);
192 Double m2(0);
193 ichi = i;
194 // Do Schwarz' calculation
195 Double b = deriv_p[iclo];
196 Double a = (deriv_p[ichi] - b) / (ichi-iclo);
197 for (Int ic=iclo; ic<=ichi; ic++) {
198 m0m += min(deriv_p[ic], 0.0);
199 Double wi = deriv_p[ic] - a*(ic-iclo) - b;
200 m0 += wi;
201 m1 += wi*ic;
202 m2 += wi*ic*ic;
203 };
204 // determinant
205 Double det = m2*m0 - m1*m1;
206 if (det > 0.0 && fabs(m0m) > FLT_EPSILON) {
207 Double xm = m1/m0;
208 Double sg = 1.69*sqrt(det) / fabs(m0);
209 // Width above critical?
210 if (sg > sigmin_p) {
211 Int is = Int(1.73*sg+0.5);
212 Int im = Int(xm+0.5);
213 Double yl(0);
214 if ((im-is) >= 0) yl = prof(im-is);
215 Double yh(0);
216 if ((im + is) <= Int(lprof_p-1)) yh = prof(im+is);
217 Double ym = prof(im);
218 // modified by dmehringer 2012apr03 to deal with 0 denominator
219 // 0.0/0.0 produces NaN on Linux but 0 on OSX
220 Double pg = (ym-0.5*(yh+yl));
221 if (pg != 0) {
222 Double denom = (1.0-exp(-0.5*(is*is)/sg/sg));
223 if (denom == 0) {
224 throw AipsError("Bailing because division by zero is undefined");
225 }
226 pg /= denom;
227 }
228 // end dmehring mods
229 pg = min(pg, ym);
230 // cout << "pg " << pg << " cutoff " << cutoff_p << endl;
231 // Above critical level? Add to list
232 if (pg > cutoff_p) {
233 // cout << pg << " " << xm << " " << sg << endl;
234 tspel.setAmpl(pg);
235 tspel.setCenter(xm);
236 tspel.setSigma(sg);
237 slist_p.insert(tspel);
238 };
239 };
240 };
241 // Next gaussian
242 iclo = ichi;
243 nmax--;
244 break;
245 }
246 default:
247 iclo = i+1;
248 break;
249 };
250 };
251}
252
253template <class MT>
254GaussianSpectralElement SpectralEstimate::convertElement (const Vector<MT>& x,
255 const GaussianSpectralElement& el) const
256{
257 GaussianSpectralElement elOut = el;
258 const Int& idxMax = x.nelements()-1;
259
260// Get current (pars are amp, center, width as the SpectralElement
261// will always be a Gaussian)
262
263 Vector<Double> par, err;
264 el.get(par);
265 el.getError(err);
266
267// Center
268
269 Int cenIdx = Int(par[1]);
270
271// Get the x-increment, local to the center, as best we can from
272// the abcissa vector. The following algorithm assumes the X
273// vector is monotonic
274
275 Double incX;
276 if (cenIdx-1<0) {
277 incX = x[1] - x[0];
278 } else if (cenIdx+1>idxMax) {
279 incX = x[idxMax] - x[idxMax-1];
280 } else {
281 incX = 0.5 * (x(cenIdx+1) - x(cenIdx-1));
282 }
283//
284 if (cenIdx<0) {
285 par[1] = incX*par[1] + x[0]; // Extrapolate from x[0]
286 } else if (cenIdx>idxMax) {
287 par[1] = incX*(par[1]-idxMax) + x[idxMax]; // Extrapolate from x[idxMax]
288 } else {
289 Double dIdx = par[1] - cenIdx;
290 par[1] = x[cenIdx] + dIdx*incX; // Interpolate
291 }
292 err[1] = abs(err[1] * incX);
293
294// Width
295
296 par[2] = abs(par[2] * incX);
297 err[2] = abs(err[2] * incX);
298
299 elOut.set(par);
300 elOut.setError(err);
301 return elOut;
302}
303
304
305} //# End namespace casa
Note: See TracBrowser for help on using the repository browser.