Changeset 1373 for trunk/src


Ignore:
Timestamp:
07/12/07 11:43:56 (17 years ago)
Author:
mar637
Message:

Added running median to smooth. This addresses Ticket #115

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/MathUtils.cpp

    r1325 r1373  
    3232#include <casa/aips.h>
    3333#include <casa/Arrays/Vector.h>
     34#include <casa/Arrays/Slice.h>
    3435#include <casa/Arrays/MaskedArray.h>
    3536#include <casa/Arrays/MaskArrMath.h>
    3637#include <casa/Arrays/VectorSTLIterator.h>
    3738#include <casa/BasicSL/String.h>
     39#include <scimath/Mathematics/MedianSlider.h>
    3840
    3941#include "MathUtils.h"
     
    147149  }
    148150}
     151
     152
     153void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
     154                             const Vector<Float>& in, const Vector<Bool>& flag,
     155                             float width)
     156{
     157  Int hwidth = Int(width+0.5);
     158  Int fwidth = hwidth*2+1;
     159  out.resize(in.nelements());
     160  outflag.resize(flag.nelements());
     161  MedianSlider ms(hwidth);
     162  Slice sl(0, fwidth-1);
     163  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
     164                  const_cast<Vector<Bool>& >(flag)(sl));
     165  uInt n = in.nelements();
     166  for (uInt i=hwidth; i<(n-hwidth); ++i) {
     167    // add data value
     168    cout << ms.nval() << endl;
     169    out[i] = ms.add(in[i], flag[i]);
     170    outflag[i] = (ms.nval() == 0);   
     171  }
     172  // replicate edge values from fisrt value with full width of values
     173  for (uInt i=0;i<hwidth;++i) {
     174    out[i] = out[hwidth];
     175    outflag[i] = outflag[hwidth];   
     176    out[n-1-i] = out[n-1-hwidth];
     177    outflag[n-1-i] = outflag[n-1-hwidth];   
     178  }
     179}
  • trunk/src/MathUtils.h

    r1325 r1373  
    5757             casa::Bool ignoreOther=casa::False);
    5858
     59/**
     60 * Apply a running median to  a masked vector.
     61 * Edge solution:  The first and last hwidth channels will be replicated
     62 * from the first/last value from a full window.
     63 * @param out the smoothed vector
     64 * @param outmask  the smoothed mask
     65 * @param in the input vector
     66 * @param mask the input mask
     67 * @param hwidth half-width of the smoothing window
     68 */
     69 void runningMedian(casa::Vector<casa::Float>& out,
     70                   casa::Vector<casa::Bool>& outmask,
     71                   const casa::Vector<casa::Float>& in,
     72                   const casa::Vector<casa::Bool>& mask,
     73                   float hwidth);
     74
    5975// Generate specified statistic
    6076float statistics(const casa::String& which,
     
    6682
    6783/**
    68  * Convert a std::vector of std::string
    69  * to a casa::Vector casa::String
    70  * @param in
    71  * @return
     84 * Convert casa implementations to stl
     85 * @param in casa string
     86 * @return a std vector of std strings
    7287 */
    7388std::vector<std::string> tovectorstring(const casa::Vector<casa::String>& in);
    7489
    7590/**
    76  * Convert a casa::Vector of casa::String
    77  * to a stl std::vector of stl std::string
     91 * convert stl implementations to casa versions
    7892 * @param in
    7993 * @return
  • trunk/src/STMath.cpp

    r1336 r1373  
    10071007}
    10081008
     1009CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
     1010                                             const std::string& kernel,
     1011                                             float width )
     1012{
     1013  CountedPtr< Scantable > out = getScantable(in, false);
     1014  Table& table = out->table();
     1015  ArrayColumn<Float> specCol(table, "SPECTRA");
     1016  ArrayColumn<uChar> flagCol(table, "FLAGTRA");
     1017  Vector<Float> spec;
     1018  Vector<uChar> flag;
     1019  for ( uInt i=0; i<table.nrow(); ++i) {
     1020    specCol.get(i, spec);
     1021    flagCol.get(i, flag);
     1022    Vector<Bool> mask(flag.nelements());
     1023    convertArray(mask, flag);
     1024    Vector<Float> specout;
     1025    Vector<Bool> maskout;
     1026    if ( kernel == "hanning" ) {
     1027      mathutil::hanning(specout, maskout, spec , !mask);
     1028      convertArray(flag, !maskout);
     1029    } else if (  kernel == "rmedian" ) {
     1030      mathutil::runningMedian(specout, maskout, spec , mask, width);
     1031      convertArray(flag, maskout);
     1032    }
     1033    flagCol.put(i, flag);
     1034    specCol.put(i, specout);
     1035  }
     1036  return out;
     1037}
     1038
    10091039CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
    10101040                                        const std::string& kernel, float width )
    10111041{
     1042  if (kernel == "rmedian"  || kernel == "hanning") {
     1043    return smoothOther(in, kernel, width);
     1044  }
    10121045  CountedPtr< Scantable > out = getScantable(in, false);
    10131046  Table& table = out->table();
     
    10321065      convertArray(mask, flag);
    10331066      Vector<Float> specout;
    1034       if ( type == VectorKernel::HANNING ) {
    1035         Vector<Bool> maskout;
    1036         mathutil::hanning(specout, maskout, spec , !mask);
    1037         convertArray(flag, !maskout);
    1038         flagCol.put(i, flag);
    1039         specCol.put(i, specout);
    1040      } else {
    1041         mathutil::replaceMaskByZero(specout, mask);
    1042         conv.linearConv(specout, spec);
    1043         specCol.put(i, specout);
    1044       }
     1067      mathutil::replaceMaskByZero(specout, mask);
     1068      conv.linearConv(specout, spec);
     1069      specCol.put(i, specout);
    10451070    }
    10461071    ++iter;
Note: See TracChangeset for help on using the changeset viewer.