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

//# Spectral2Estimate.cc: Member templates for SpectralEstimate
//# Copyright (C) 2001,2002,2003,2004
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//# $Id: Spectral2Estimate.tcc 21465 2014-06-19 05:56:56Z gervandiepen $

//# Includes
#include <components/SpectralComponents/SpectralEstimate.h>

#include <casa/BasicMath/Math.h>
#include <casa/BasicSL/Constants.h>
#include <casa/Utilities/Assert.h>
#include <components/SpectralComponents/CompiledSpectralElement.h>
#include <components/SpectralComponents/GaussianSpectralElement.h>
#include <components/SpectralComponents/PolynomialSpectralElement.h>

namespace casa { //#Begin namespace casa

//# Member templates
template <class MT>
const SpectralList &SpectralEstimate::estimate(const Vector<MT> &prof,
                                               Vector<MT> *der) {
          if (prof.nelements() != lprof_p) {
    delete [] deriv_p; deriv_p = 0; lprof_p = 0;
    lprof_p = prof.nelements();
    deriv_p = new Double[lprof_p];
  };
  // Check if signal in window
  if (!window(prof)) return slist_p;
  // Limit window
  windowEnd_p = min(windowEnd_p+q_p , Int(lprof_p)); 
  windowLow_p = max(windowLow_p-q_p , 0 );
  // Get the second derivatives
  findc2(prof);
  // Next for debugging
  if (der) {
    for (uInt i=0; i<lprof_p; i++) (*der)[i] = deriv_p[i];
  };
  // Find the estimates (sorted)
  findga(prof);
  // cout << slist_p << endl;
  return slist_p;
}

template <class MT>
const SpectralList& SpectralEstimate::estimate(const Vector<MT>& x,
                                               const Vector<MT>& y)
{
  if (x.nelements() != y.nelements()) {
     throw(AipsError("Abcissa and ordinate vectors must be the same length"));
  }
  if (x.nelements()==1) {
     throw(AipsError("Not enough elements in vectors")); 
  }
  // Get pixel-based estimate (into slist_p)
  estimate(y);
  // Convert
  for (uInt i=0; i<slist_p.nelements(); i++) {
          if (slist_p[i]->getType() != SpectralElement::GAUSSIAN) {
                  throw AipsError("Non-gaussian spectral types cannot be estimated");
          }
          const GaussianSpectralElement elIn = *dynamic_cast<const GaussianSpectralElement *>(slist_p[i]);
          GaussianSpectralElement elOut = convertElement (x, elIn);
          slist_p.set(elOut, i);
  }
  return slist_p;
}


template <class MT>
uInt SpectralEstimate::window(const Vector<MT> &prof) {
  windowLow_p =0;
  windowEnd_p = 0;
  if (!useWindow_p || rms_p <= 0.0 || lprof_p == 0) {
    if (regionEnd_p) {
      windowLow_p = min(max(0,regionLow_p),Int(lprof_p));
      windowEnd_p = min(regionEnd_p, Int(lprof_p));
    } else windowEnd_p = lprof_p;
    return windowEnd_p-windowLow_p;
  };
  // Total flux in profile and max position
  Double flux(0.0);
  Double pmax(prof(0));
  uInt imax(0);
  for (Int i=windowLow_p; i<windowEnd_p; i++) {
    if (prof(i)>pmax) {
      pmax = prof(i);
      imax = i;
    };
    flux += prof(i);
  };
  // No data
  if (pmax < cutoff_p) return 0;
  // Window boundaries; new/old base and centre; width
  Int width(-1);
  Int nw(0);
  Double bnew(flux), bold;
  Double cnew(imax), cold;
  do {
    width++;
    cold = cnew;
    bold = bnew;
    windowLow_p = max(0, Int(cold-width+0.5));
    windowEnd_p = min(Int(lprof_p), Int(cold+width+1.5));
    // flux and first moment in window
    Double s(0);
    Double c(0);
    for (Int i=windowLow_p; i<windowEnd_p; i++) {
      s += prof(i);
      c += i*prof(i);
    };
    bnew = flux-s;
    nw = lprof_p-windowEnd_p+windowLow_p;
    if (s != 0.0) {
      cnew = c/s;
      if (cnew < 0 || cnew >= lprof_p) cnew = cold;
    };
  } while (abs(bnew-bold) > rms_p && nw);
  return windowEnd_p-windowLow_p;
}

template <class MT>
void SpectralEstimate::findc2(const Vector<MT> &prof) {
  for (Int i=windowLow_p; i<windowEnd_p; i++) {
    // Moments
    Double m0(0.0); 
    Double m2(0.0); 
    for (Int j = -q_p; j <= q_p; j++) {
      Int k = i+j;
      if (k >= 0 && k<Int(lprof_p)) {
        // add to moments
        m0 += prof(k);
        m2 += prof(k)*j*j;
      };
    };
    // get the derivative
    deriv_p[i] = a_p*(m2-b_p*m0);
  };
}

template <class MT>
void SpectralEstimate::findga(const Vector<MT> &prof) {
        Int i(windowLow_p-1);
        // Window on Gaussian
        Int iclo(windowLow_p);
        Int ichi(windowLow_p);
        // Peak counter
        Int nmax = 0;
        GaussianSpectralElement tspel;
        while (++i < windowEnd_p) {
        if (deriv_p[i] > 0.0) {
                        // At edge?
                        if (i > windowLow_p && i < windowEnd_p-1) {
                                // Peak in 2nd derivative
                                if (deriv_p[i-1] < deriv_p[i] && deriv_p[i+1] < deriv_p[i]) nmax++;
                                // At start
                        } else if (i == windowLow_p && deriv_p[i+1] < deriv_p[i]) nmax++;
                        // At end of window
                        else if (i == windowEnd_p-1 && deriv_p[i-1] < deriv_p[i]) nmax++;
                };
                switch (nmax) {
                // Search for next peak
                case 1:
                        break;
                        // Found a Gaussian
                case 2: {
                        // Some moments
                        Double m0m(0);
                        Double m0(0);
                        Double m1(0);
                        Double m2(0);
                        ichi = i;
                        // Do Schwarz' calculation
                        Double b = deriv_p[iclo];
                        Double a = (deriv_p[ichi] - b) / (ichi-iclo);
                        for (Int ic=iclo; ic<=ichi; ic++) {
                                m0m += min(deriv_p[ic], 0.0);
                                Double wi = deriv_p[ic] - a*(ic-iclo) - b;
                                m0 += wi;
                                m1 += wi*ic;
                                m2 += wi*ic*ic;
                        };
                        // determinant
                        Double det = m2*m0 - m1*m1;
                        if (det > 0.0 && fabs(m0m) >  FLT_EPSILON) {
                                Double   xm = m1/m0;
                                Double sg = 1.69*sqrt(det) / fabs(m0);
                                // Width above critical?
                                if (sg > sigmin_p) {
                                        Int is = Int(1.73*sg+0.5);
                                        Int im = Int(xm+0.5);
                                        Double yl(0);
                                        if ((im-is) >= 0) yl = prof(im-is);
                                        Double yh(0);
                                        if ((im + is) <= Int(lprof_p-1)) yh = prof(im+is);
                                        Double ym = prof(im);
                    // modified by dmehringer 2012apr03 to deal with 0 denominator
                    // 0.0/0.0 produces NaN on Linux but 0 on OSX
                    Double pg = (ym-0.5*(yh+yl));
                                        if (pg != 0) {
                                                Double denom = (1.0-exp(-0.5*(is*is)/sg/sg));
                                                if (denom == 0) {
                                                        throw AipsError("Bailing because division by zero is undefined");
                                                }
                                                pg /= denom;
                                        }
                    // end dmehring mods
                                        pg = min(pg, ym);
                                        // cout << "pg " << pg << " cutoff " << cutoff_p << endl;
                                        // Above critical level? Add to list
                                        if (pg > cutoff_p) {
                                        //      cout << pg << " " << xm << " " << sg << endl;
                                                tspel.setAmpl(pg);
                                                tspel.setCenter(xm);
                                                tspel.setSigma(sg);
                        slist_p.insert(tspel);
                                        };
                                };
                        };
                        // Next gaussian
                        iclo = ichi;
                        nmax--;
                        break;
                }
                default:
                        iclo = i+1;
                        break;
                };
        };
}

template <class MT>
GaussianSpectralElement SpectralEstimate::convertElement (const Vector<MT>& x,
                                                  const GaussianSpectralElement& el) const
{
        GaussianSpectralElement elOut = el;
   const Int& idxMax = x.nelements()-1;

// Get current (pars are amp, center, width as the SpectralElement
// will always be a Gaussian)

   Vector<Double> par, err;
   el.get(par);
   el.getError(err);

// Center

   Int cenIdx = Int(par[1]);

// Get the x-increment, local to the center, as best we can from 
// the abcissa vector.  The following algorithm assumes the X
// vector is monotonic

   Double incX;
   if (cenIdx-1<0) {
      incX = x[1] - x[0];
   } else if (cenIdx+1>idxMax) {
      incX = x[idxMax] - x[idxMax-1];
   } else {
      incX = 0.5 * (x(cenIdx+1) - x(cenIdx-1));  
   }
//
   if (cenIdx<0) {
      par[1] = incX*par[1] + x[0];                 // Extrapolate from x[0]
   } else if (cenIdx>idxMax) {
      par[1] = incX*(par[1]-idxMax) + x[idxMax];   // Extrapolate from x[idxMax]
   } else {
      Double dIdx = par[1] - cenIdx;
      par[1] = x[cenIdx] + dIdx*incX;              // Interpolate
   }
   err[1] = abs(err[1] * incX);

// Width

   par[2] = abs(par[2] * incX);
   err[2] = abs(err[2] * incX);

   elOut.set(par);
   elOut.setError(err);
   return elOut;
}


} //# End namespace casa