Line data    Source code

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