Changeset 1373

Timestamp:

07/12/07 11:43:56 (17 years ago)

Author:

mar637

Message:

Added running median to smooth. This addresses Ticket #115

Location:

trunk

Files:

• python/scantable.py (modified) (1 diff)
• src/MathUtils.cpp (modified) (2 diffs)
• src/MathUtils.h (modified) (2 diffs)
• src/STMath.cpp (modified) (2 diffs)

trunk/python/scantable.py

@@ -1192 +1192 @@
         Parameters:
             kernel:     The type of smoothing kernel. Select from
-                        'hanning' (default), 'gaussian' and 'boxcar'.
-                        The first three characters are sufficient.
+                        'hanning' (default), 'gaussian', 'boxcar' and
+                        'rmedian'
             width:      The width of the kernel in pixels. For hanning this is
                         ignored otherwise it defauls to 5 pixels.
                         For 'gaussian' it is the Full Width Half
                         Maximum. For 'boxcar' it is the full width.
+                        For 'rmedian' it is the half width.
             insitu:     if False a new scantable is returned.
                         Otherwise, the scaling is done in-situ

trunk/src/MathUtils.cpp

@@ -32 +32 @@
 #include <casa/aips.h>
 #include <casa/Arrays/Vector.h>
+#include <casa/Arrays/Slice.h>
 #include <casa/Arrays/MaskedArray.h>
 #include <casa/Arrays/MaskArrMath.h>
 #include <casa/Arrays/VectorSTLIterator.h>
 #include <casa/BasicSL/String.h>
+#include <scimath/Mathematics/MedianSlider.h>
 
 #include "MathUtils.h"
@@ -147 +149 @@
   }
 }
+
+
+void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
+                             const Vector<Float>& in, const Vector<Bool>& flag,
+                             float width)
+{
+  Int hwidth = Int(width+0.5);
+  Int fwidth = hwidth*2+1;
+  out.resize(in.nelements());
+  outflag.resize(flag.nelements());
+  MedianSlider ms(hwidth);
+  Slice sl(0, fwidth-1);
+  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
+    cout << ms.nval() << endl;
+    out[i] = ms.add(in[i], flag[i]); 
+    outflag[i] = (ms.nval() == 0);    
+  }
+  // replicate edge values from fisrt value with full width of values
+  for (uInt i=0;i<hwidth;++i) {
+    out[i] = out[hwidth];
+    outflag[i] = outflag[hwidth];    
+    out[n-1-i] = out[n-1-hwidth];
+    outflag[n-1-i] = outflag[n-1-hwidth];    
+  }
+}

trunk/src/MathUtils.h

@@ -57 +57 @@
              casa::Bool ignoreOther=casa::False);
 
+/**
+ * Apply a running median to  a masked vector.
+ * 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
+ * @param outmask  the smoothed mask
+ * @param in the input vector
+ * @param mask the input mask
+ * @param hwidth half-width of the smoothing window
+ */
+ void runningMedian(casa::Vector<casa::Float>& out, 
+                   casa::Vector<casa::Bool>& outmask,
+                   const casa::Vector<casa::Float>& in, 
+                   const casa::Vector<casa::Bool>& mask,
+                   float hwidth);
+
 // Generate specified statistic
 float statistics(const casa::String& which,
@@ -66 +82 @@
 
 /**
- * Convert a std::vector of std::string
- * to a casa::Vector casa::String
- * @param in
- * @return
+ * Convert casa implementations to stl
+ * @param in casa string
+ * @return a std vector of std strings
  */
 std::vector<std::string> tovectorstring(const casa::Vector<casa::String>& in);
 
 /**
- * Convert a casa::Vector of casa::String
- * to a stl std::vector of stl std::string
+ * convert stl implementations to casa versions
  * @param in
  * @return

trunk/src/STMath.cpp

@@ -1007 +1007 @@
 }
 
+CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
+                                             const std::string& kernel,
+                                             float width )
+{
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& table = out->table();
+  ArrayColumn<Float> specCol(table, "SPECTRA");
+  ArrayColumn<uChar> flagCol(table, "FLAGTRA");
+  Vector<Float> spec;
+  Vector<uChar> flag;
+  for ( uInt i=0; i<table.nrow(); ++i) {
+    specCol.get(i, spec);
+    flagCol.get(i, flag);
+    Vector<Bool> mask(flag.nelements());
+    convertArray(mask, flag);
+    Vector<Float> specout;
+    Vector<Bool> maskout;
+    if ( kernel == "hanning" ) {
+      mathutil::hanning(specout, maskout, spec , !mask);
+      convertArray(flag, !maskout);
+    } else if (  kernel == "rmedian" ) {
+      mathutil::runningMedian(specout, maskout, spec , mask, width);
+      convertArray(flag, maskout);
+    }
+    flagCol.put(i, flag);
+    specCol.put(i, specout);
+  }
+  return out;
+}
+
 CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
                                         const std::string& kernel, float width )
 {
+  if (kernel == "rmedian"  || kernel == "hanning") {
+    return smoothOther(in, kernel, width);
+  }
   CountedPtr< Scantable > out = getScantable(in, false);
   Table& table = out->table();
@@ -1032 +1065 @@
       convertArray(mask, flag);
       Vector<Float> specout;
-      if ( type == VectorKernel::HANNING ) {
-        Vector<Bool> maskout;
-        mathutil::hanning(specout, maskout, spec , !mask);
-        convertArray(flag, !maskout);
-        flagCol.put(i, flag);
-        specCol.put(i, specout);
-     } else {
-        mathutil::replaceMaskByZero(specout, mask);
-        conv.linearConv(specout, spec);
-        specCol.put(i, specout);
-      }
+      mathutil::replaceMaskByZero(specout, mask);
+      conv.linearConv(specout, spec);
+      specCol.put(i, specout);
     }
     ++iter;