[2980] | 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 | |
---|
| 38 | namespace casa { //#Begin namespace casa |
---|
| 39 | |
---|
| 40 | //# Member templates |
---|
| 41 | template <class MT> |
---|
| 42 | const 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 | |
---|
| 66 | template <class MT> |
---|
| 67 | const 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 | |
---|
| 91 | template <class MT> |
---|
| 92 | uInt 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 | |
---|
| 143 | template <class MT> |
---|
| 144 | void 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 | |
---|
| 162 | template <class MT> |
---|
| 163 | void 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 | |
---|
| 253 | template <class MT> |
---|
| 254 | GaussianSpectralElement 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 |
---|