Changeset 1570 for trunk


Ignore:
Timestamp:
06/29/09 12:04:00 (15 years ago)
Author:
Malte Marquarding
Message:

Ticket #167: c++ part of running polynomial smoothing

Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/MathUtils.cpp

    r1412 r1570  
    3838#include <casa/BasicSL/String.h>
    3939#include <scimath/Mathematics/MedianSlider.h>
     40
     41#include <scimath/Fitting/LinearFit.h>
     42#include <scimath/Functionals/Polynomial.h>
     43#include <scimath/Mathematics/AutoDiff.h>
     44
    4045
    4146#include "MathUtils.h"
     
    161166  MedianSlider ms(hwidth);
    162167  Slice sl(0, fwidth-1);
    163   Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl), 
     168  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
    164169                  const_cast<Vector<Bool>& >(flag)(sl));
    165170  uInt n = in.nelements();
    166171  for (uInt i=hwidth; i<(n-hwidth); ++i) {
    167172    // add data value
    168     out[i] = ms.add(in[i+hwidth], flag[i+hwidth]); 
    169     outflag[i] = (ms.nval() == 0);   
    170   }
    171   // replicate edge values from fisrt value with full width of values
     173    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
     174    outflag[i] = (ms.nval() == 0);
     175  }
     176  // replicate edge values from first value with full width of values
    172177  for (uInt i=0;i<hwidth;++i) {
    173178    out[i] = out[hwidth];
    174     outflag[i] = outflag[hwidth];   
     179    outflag[i] = outflag[hwidth];
    175180    out[n-1-i] = out[n-1-hwidth];
    176     outflag[n-1-i] = outflag[n-1-hwidth];   
    177   }
    178 }
     181    outflag[n-1-i] = outflag[n-1-hwidth];
     182  }
     183}
     184
     185void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
     186                       const Vector<Float>& in, const Vector<Bool>& mask,
     187                       float width, int order)
     188{
     189  Int hwidth = Int(width+0.5);
     190  Int fwidth = hwidth*2+1;
     191  out.resize(in.nelements());
     192  outmask.resize(mask.nelements());
     193  LinearFit<Float> fitter;
     194  Polynomial<Float> poly(order);
     195  fitter.setFunction(poly);
     196  Vector<Float> sigma(fwidth);
     197  sigma = 1.0;
     198  Vector<Float> parms;
     199  Vector<Float> x(fwidth);
     200  indgen(x);
     201
     202  uInt n = in.nelements();
     203
     204  for (uInt i=hwidth; i<(n-hwidth); ++i) {
     205    // add data value
     206    if (mask[i]) {
     207      Slice sl(i-hwidth, fwidth);
     208      const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
     209      const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
     210      parms = fitter.fit(x, y, sigma, &m);
     211
     212      poly.setCoefficients(parms);
     213      out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
     214    } else {
     215      out[i] = in[i];
     216    }
     217    outmask[i] = mask[i];
     218  }
     219  // replicate edge values from first value with full width of values
     220  for (uInt i=0;i<hwidth;++i) {
     221    out[i] = out[hwidth];
     222    outmask[i] = outmask[hwidth];
     223    out[n-1-i] = out[n-1-hwidth];
     224    outmask[n-1-i] = outmask[n-1-hwidth];
     225  }
     226}
  • trunk/src/MathUtils.h

    r1373 r1570  
    5050 * @param ignoreOther drop every second channel (NYI)
    5151 */
    52 void hanning(casa::Vector<casa::Float>& out, 
     52void hanning(casa::Vector<casa::Float>& out,
    5353             casa::Vector<casa::Bool>& outmask,
    54              const casa::Vector<casa::Float>& in, 
     54             const casa::Vector<casa::Float>& in,
    5555             const casa::Vector<casa::Bool>& mask,
    5656             casa::Bool relaxed=casa::False,
     
    5959/**
    6060 * Apply a running median to  a masked vector.
    61  * Edge solution:  The first and last hwidth channels will be replicated 
     61 * Edge solution:  The first and last hwidth channels will be replicated
    6262 * from the first/last value from a full window.
    6363 * @param out the smoothed vector
     
    6767 * @param hwidth half-width of the smoothing window
    6868 */
    69  void runningMedian(casa::Vector<casa::Float>& out,
     69 void runningMedian(casa::Vector<casa::Float>& out,
     70                   casa::Vector<casa::Bool>& outflag,
     71                   const casa::Vector<casa::Float>& in,
     72                   const casa::Vector<casa::Bool>& flag,
     73                   float hwidth);
     74
     75 void polyfit(casa::Vector<casa::Float>& out,
    7076                   casa::Vector<casa::Bool>& outmask,
    71                    const casa::Vector<casa::Float>& in, 
     77                   const casa::Vector<casa::Float>& in,
    7278                   const casa::Vector<casa::Bool>& mask,
    73                    float hwidth);
     79                   float hwidth, int order);
     80
    7481
    7582// Generate specified statistic
  • trunk/src/STMath.cpp

    r1569 r1570  
    13681368    if (d < 0) {
    13691369      Instrument inst =
    1370         STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 
     1370        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
    13711371                                  True);
    13721372      STAttr sda;
     
    14571457CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
    14581458                                             const std::string& kernel,
    1459                                              float width )
     1459                                             float width, int order)
    14601460{
    14611461  CountedPtr< Scantable > out = getScantable(in, false);
     
    14781478      mathutil::runningMedian(specout, maskout, spec , mask, width);
    14791479      convertArray(flag, maskout);
     1480    } else if ( kernel == "poly" ) {
     1481      mathutil::polyfit(specout, maskout, spec, !mask, width, order);
     1482      convertArray(flag, !maskout);
    14801483    }
    14811484    flagCol.put(i, flag);
  • trunk/src/STMath.h

    r1391 r1570  
    120120
    121121  casa::CountedPtr<Scantable>
    122     binaryOperate( const casa::CountedPtr<Scantable>& left, 
    123                    const casa::CountedPtr<Scantable>& right, 
     122    binaryOperate( const casa::CountedPtr<Scantable>& left,
     123                   const casa::CountedPtr<Scantable>& right,
    124124                   const std::string& mode);
    125125
     
    215215  casa::CountedPtr<Scantable>
    216216    smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel,
    217                       float width);
     217                      float width, int order=2);
    218218
    219219  casa::CountedPtr<Scantable>
     
    263263   */
    264264  casa::CountedPtr<Scantable>
    265     lagFlag( const casa::CountedPtr<Scantable>& in, double frequency,
    266               double width);
     265    lagFlag( const casa::CountedPtr<Scantable>& in, double start,
     266             double end, const std::string& mode="frequency");
    267267
    268268private:
     
    290290                              bool tokelvin, float cfac);
    291291
    292   casa::CountedPtr< Scantable > 
     292  casa::CountedPtr< Scantable >
    293293    smoothOther( const casa::CountedPtr< Scantable >& in,
    294294                 const std::string& kernel,
    295                  float width );
     295                 float width, int order=2 );
    296296
    297297  casa::CountedPtr< Scantable >
  • trunk/src/STMathWrapper.h

    r1541 r1570  
    130130
    131131  ScantableWrapper
    132     smooth(const ScantableWrapper& in, const std::string& kernel, float width)
    133   { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width)); }
     132    smooth(const ScantableWrapper& in, const std::string& kernel, float width,
     133           int order=2)
     134  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width, order)); }
    134135
    135136  ScantableWrapper
     
    185186
    186187  ScantableWrapper lagFlag( const ScantableWrapper& in,
    187                             double frequency, double width )
    188   { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
     188                            double start, double end,
     189                            const std::string& mode="frequency" )
     190  { return ScantableWrapper(STMath::lagFlag(in.getCP(), start, end,
     191                                            mode)); }
    189192
    190193};
Note: See TracChangeset for help on using the changeset viewer.