Changeset 1570

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:

• MathUtils.cpp (modified) (2 diffs)
• MathUtils.h (modified) (3 diffs)
• STMath.cpp (modified) (3 diffs)
• STMath.h (modified) (4 diffs)
• STMathWrapper.h (modified) (2 diffs)

trunk/src/MathUtils.cpp

@@ -38 +38 @@
 #include <casa/BasicSL/String.h>
 #include <scimath/Mathematics/MedianSlider.h>
+
+#include <scimath/Fitting/LinearFit.h>
+#include <scimath/Functionals/Polynomial.h>
+#include <scimath/Mathematics/AutoDiff.h>
+
 
 #include "MathUtils.h"
@@ -161 +166 @@
   MedianSlider ms(hwidth);
   Slice sl(0, fwidth-1);
-  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl), 
+  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
                   const_cast<Vector<Bool>& >(flag)(sl));
   uInt n = in.nelements();
   for (uInt i=hwidth; i<(n-hwidth); ++i) {
     // add data value
-    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]); 
-    outflag[i] = (ms.nval() == 0);    
-  }
-  // replicate edge values from fisrt value with full width of values
+    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
+    outflag[i] = (ms.nval() == 0);
+  }
+  // replicate edge values from first value with full width of values
   for (uInt i=0;i<hwidth;++i) {
     out[i] = out[hwidth];
-    outflag[i] = outflag[hwidth];    
+    outflag[i] = outflag[hwidth];
     out[n-1-i] = out[n-1-hwidth];
-    outflag[n-1-i] = outflag[n-1-hwidth];    
-  }
-}
+    outflag[n-1-i] = outflag[n-1-hwidth];
+  }
+}
+
+void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
+                       const Vector<Float>& in, const Vector<Bool>& mask,
+                       float width, int order)
+{
+  Int hwidth = Int(width+0.5);
+  Int fwidth = hwidth*2+1;
+  out.resize(in.nelements());
+  outmask.resize(mask.nelements());
+  LinearFit<Float> fitter;
+  Polynomial<Float> poly(order);
+  fitter.setFunction(poly);
+  Vector<Float> sigma(fwidth);
+  sigma = 1.0;
+  Vector<Float> parms;
+  Vector<Float> x(fwidth);
+  indgen(x);
+
+  uInt n = in.nelements();
+
+  for (uInt i=hwidth; i<(n-hwidth); ++i) {
+    // add data value
+    if (mask[i]) {
+      Slice sl(i-hwidth, fwidth);
+      const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
+      const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
+      parms = fitter.fit(x, y, sigma, &m);
+
+      poly.setCoefficients(parms);
+      out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
+    } else {
+      out[i] = in[i];
+    }
+    outmask[i] = mask[i];
+  }
+  // replicate edge values from first value with full width of values
+  for (uInt i=0;i<hwidth;++i) {
+    out[i] = out[hwidth];
+    outmask[i] = outmask[hwidth];
+    out[n-1-i] = out[n-1-hwidth];
+    outmask[n-1-i] = outmask[n-1-hwidth];
+  }
+}

trunk/src/MathUtils.h

@@ -50 +50 @@
  * @param ignoreOther drop every second channel (NYI)
  */
-void hanning(casa::Vector<casa::Float>& out, 
+void hanning(casa::Vector<casa::Float>& out,
              casa::Vector<casa::Bool>& outmask,
-             const casa::Vector<casa::Float>& in, 
+             const casa::Vector<casa::Float>& in,
              const casa::Vector<casa::Bool>& mask,
              casa::Bool relaxed=casa::False,
@@ -59 +59 @@
 /**
  * Apply a running median to  a masked vector.
- * Edge solution:  The first and last hwidth channels will be replicated 
+ * Edge solution:  The first and last hwidth channels will be replicated
  * from the first/last value from a full window.
  * @param out the smoothed vector
@@ -67 +67 @@
  * @param hwidth half-width of the smoothing window
  */
- void runningMedian(casa::Vector<casa::Float>& out, 
+ void runningMedian(casa::Vector<casa::Float>& out,
+                   casa::Vector<casa::Bool>& outflag,
+                   const casa::Vector<casa::Float>& in,
+                   const casa::Vector<casa::Bool>& flag,
+                   float hwidth);
+
+ void polyfit(casa::Vector<casa::Float>& out,
                    casa::Vector<casa::Bool>& outmask,
-                   const casa::Vector<casa::Float>& in, 
+                   const casa::Vector<casa::Float>& in,
                    const casa::Vector<casa::Bool>& mask,
-                   float hwidth);
+                   float hwidth, int order);
+
 
 // Generate specified statistic

trunk/src/STMath.cpp

@@ -1368 +1368 @@
     if (d < 0) {
       Instrument inst =
-        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 
+        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
                                   True);
       STAttr sda;
@@ -1457 +1457 @@
 CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
                                              const std::string& kernel,
-                                             float width )
+                                             float width, int order)
 {
   CountedPtr< Scantable > out = getScantable(in, false);
@@ -1478 +1478 @@
       mathutil::runningMedian(specout, maskout, spec , mask, width);
       convertArray(flag, maskout);
+    } else if ( kernel == "poly" ) {
+      mathutil::polyfit(specout, maskout, spec, !mask, width, order);
+      convertArray(flag, !maskout);
     }
     flagCol.put(i, flag);

trunk/src/STMath.h

@@ -120 +120 @@
 
   casa::CountedPtr<Scantable>
-    binaryOperate( const casa::CountedPtr<Scantable>& left, 
-                   const casa::CountedPtr<Scantable>& right, 
+    binaryOperate( const casa::CountedPtr<Scantable>& left,
+                   const casa::CountedPtr<Scantable>& right,
                    const std::string& mode);
 
@@ -215 +215 @@
   casa::CountedPtr<Scantable>
     smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel,
-                      float width);
+                      float width, int order=2);
 
   casa::CountedPtr<Scantable>
@@ -263 +263 @@
    */
   casa::CountedPtr<Scantable>
-    lagFlag( const casa::CountedPtr<Scantable>& in, double frequency,
-              double width);
+    lagFlag( const casa::CountedPtr<Scantable>& in, double start,
+             double end, const std::string& mode="frequency");
 
 private:
@@ -290 +290 @@
                               bool tokelvin, float cfac);
 
-  casa::CountedPtr< Scantable > 
+  casa::CountedPtr< Scantable >
     smoothOther( const casa::CountedPtr< Scantable >& in,
                  const std::string& kernel,
-                 float width );
+                 float width, int order=2 );
 
   casa::CountedPtr< Scantable >

trunk/src/STMathWrapper.h

@@ -130 +130 @@
 
   ScantableWrapper
-    smooth(const ScantableWrapper& in, const std::string& kernel, float width)
-  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width)); }
+    smooth(const ScantableWrapper& in, const std::string& kernel, float width,
+           int order=2)
+  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width, order)); }
 
   ScantableWrapper
@@ -185 +186 @@
 
   ScantableWrapper lagFlag( const ScantableWrapper& in,
-                            double frequency, double width )
-  { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
+                            double start, double end,
+                            const std::string& mode="frequency" )
+  { return ScantableWrapper(STMath::lagFlag(in.getCP(), start, end,
+                                            mode)); }
 
 };