Changeset 2125


Ignore:
Timestamp:
04/08/11 21:01:15 (14 years ago)
Author:
WataruKawasaki
Message:

New Development: No

JIRA Issue: Yes CAS-2776

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): SD

Description: bugfix of STMath::smoothOther(), STMath::average(), and RowAccumulator (used in sdsmooth)


Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/MathUtils.cpp

    r1819 r2125  
    3030//#---------------------------------------------------------------------------
    3131
     32#include <iostream>
    3233#include <casa/aips.h>
    3334#include <casa/Arrays/Vector.h>
    … …  
    161162  out.resize(in.nelements());
    162163  out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
     164  //--s
     165  cout << "out[0]=" << out[0] << ";" << endl;
     166  //--e
    163167  outmask.resize(mask.nelements());
    164   outmask = False;
     168  outmask[0] = mask[0]; outmask[outmask.nelements()-1] = mask[mask.nelements()-1];
    165169  uInt m;Vector<Float>* w;
    166170  for (uInt i=1; i < out.nelements()-1;++i) {
    … …  
    169173    if (weighted[m]) {
    170174      out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
    171       outmask[i] = True;
     175      //outmask[i] = True;
    172176    } else { // mask it
    173177      out[i] = in[i];//use arbitrary value
    174       outmask[i] = False;
     178      //outmask[i] = False;
    175179    }
     180    //--s
     181    if (i < 13) {
     182      cout << "mask[" << i-1 << "-" << i+1 << "]=[" << mask[i-1] << ", " << mask[i] << ", " << mask[i+1] << "];  ";
     183      cout << "m="   << m   << ";  ";
     184      cout << "weights[" << m << "]="   << weights[m]   << ";  ";
     185      cout << "in[" << i-1 << "-" << i+1 << "]=[" << in[i-1] << ", " << in[i] << ", " << in[i+1] << "] --> ";
     186      cout << out[i]  << ";" << endl;
     187      cout << "-------------------------------" << endl;
     188    }
     189    //--e
     190    outmask[i] = mask[i];
    176191  }
    177192}

  • trunk/src/RowAccumulator.cpp

    r1819 r2125  
    1010//
    1111//
     12
     13#include <iostream>
     14
    1215#include <casa/iomanip.h>
    13 
    1416#include <casa/Arrays/MaskArrMath.h>
    1517#include <casa/Arrays/MaskArrLogi.h>
    … …  
    1820#include "RowAccumulator.h"
    1921
    20 
    2122using namespace casa;
    2223using namespace asap;
    2324
    24 RowAccumulator::RowAccumulator(WeightType wt) :
    25   weightType_(wt),
    26   initialized_(False)
     25RowAccumulator::RowAccumulator(WeightType wt) : weightType_(wt), initialized_(False)
    2726{
    2827  reset();
    … …  
    3433
    3534
    36 void RowAccumulator::add( const Vector< Float >& v,
    37                           const Vector< Bool >& m,
    38                           const Vector< Float >& tsys,
    39                           Double interval,
    40                           Double time )
    41 {
    42   if (!initialized_) {
    43     Vector<Float> dummy(v.nelements(), 0.0);
    44     Vector<Bool> dummymsk(m.nelements(), True);
    45     spectrum_.setData(dummy, dummymsk);
    46     n_.setData(Vector<uInt>(v.nelements(), 0), dummymsk);
    47     weightSum_.setData(Vector<Float>(v.nelements(), 0.0), dummymsk);
    48     tsysSum_.resize(tsys.nelements()); tsysSum_=0.0;
    49   }
    50   // add spectrum related weights, so far it is variance only.
    51   Float totalweight = 1.0;
    52 
    53   // only add these if not everything masked
    54   if ( !allEQ(m, False) ) {
    55     totalweight *= addTsys(tsys);
    56     totalweight *= addInterval(interval);
    57     addTime(time);
    58   }
    59   addSpectrum(v, m, totalweight);
    60   initialized_ = True;
    61 }
    62 
    63 void RowAccumulator::addSpectrum( const Vector< Float >& v,
    64                                   const Vector< Bool >& m,
    65                                   Float weight)
    66 {
    67   Float totalweight = weight;
    68   MaskedArray<Float> data(v,m);
    69   if ( weightType_ == asap::W_VAR ) {
    70     if (m.nelements() == userMask_.nelements()) {
    71       Float fac = 1.0/variance(data(userMask_));
    72       totalweight *= fac;
    73     } else {
    74       Float fac = 1.0/variance(data);
    75       totalweight *= fac;
    76     }
    77   }
    78   data *= totalweight;
    79   MaskedArray<Float> wadd(Vector<Float>(m.nelements(),totalweight), m);
    80   weightSum_ += wadd;
    81   spectrum_ += data;
    82   const MaskedArray<uInt> inc(Vector<uInt>(m.nelements(),1), m);
    83   n_ += inc;
    84 }
    85 
    86 Float RowAccumulator::addTsys( const casa::Vector< casa::Float > & v )
     35void RowAccumulator::reset(const uInt size, const uInt tsysSize)
     36{
     37  Vector<Bool> m(size, True);
     38
     39  spectrum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
     40  spectrumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
     41
     42  n_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
     43  nNoMask_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
     44
     45  weightSum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
     46  weightSumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
     47
     48  tsysSum_.resize(tsysSize); tsysSum_=0.0;
     49  tsysSumNoMask_.resize(tsysSize); tsysSumNoMask_=0.0;
     50
     51  intervalSum_ = 0.0;
     52  intervalSumNoMask_ = 0.0;
     53
     54  timeSum_ = 0.0;
     55  timeSumNoMask_ = 0.0;
     56}
     57
     58void RowAccumulator::initialize(const uInt size, const uInt tsysSize)
     59{
     60    reset(size, tsysSize);
     61    initialized_ = True;
     62}
     63
     64void RowAccumulator::add(const Vector<Float>& v,
     65                         const Vector<Bool>& m,
     66                         const Vector<Float>& tsys,
     67                         const Double interval,
     68                         const Double time)
     69{
     70  uInt size = v.nelements();
     71  //if (size != m.nelements()) raiseError;
     72  if (!initialized_) initialize(size, tsys.nelements());
     73
     74  addSpectrum(v, m, tsys, interval, time);
     75}
     76
     77void RowAccumulator::addSpectrum(const Vector<Float>& v,
     78                                 const Vector<Bool>& m,
     79                                 const Vector<Float>& tsys,
     80                                 const Double interval,
     81                                 const Double time)
     82{
     83  doAddSpectrum(v, m, tsys, interval, time, False);
     84  doAddSpectrum(v, m, tsys, interval, time, True);  // CAS-2776
     85}
     86
     87void RowAccumulator::doAddSpectrum(const Vector<Float>& v,
     88                                   const Vector<Bool>& m,
     89                                   const Vector<Float>& tsys,
     90                                   const Double interval,
     91                                   const Double time,
     92                                   const Bool ignoreMask)
     93{
     94  Vector<Float> vUse = v.copy();
     95  Vector<Bool> mUse = m.copy();
     96  if (ignoreMask) mUse = !mUse;
     97
     98  MaskedArray<Float> vadd(vUse, mUse);
     99  Float totalWeight = getTotalWeight(vadd, tsys, interval, time, ignoreMask);
     100  vadd *= totalWeight;
     101  MaskedArray<Float> wadd(Vector<Float>(mUse.nelements(), totalWeight), mUse);
     102  MaskedArray<uInt> inc(Vector<uInt>(mUse.nelements(), 1), mUse);
     103
     104  if (ignoreMask) {
     105    spectrumNoMask_ += vadd;
     106    weightSumNoMask_ += wadd;
     107    nNoMask_ += inc;
     108  } else {
     109    spectrum_ += vadd;
     110    weightSum_ += wadd;
     111    n_ += inc;
     112  }
     113
     114  //
     115  cout << "***" << endl;
     116  Vector<Float> spe = spectrum_.getArray();
     117  Vector<Float> spe0 = spectrumNoMask_.getArray();
     118  Vector<Float> wei = weightSum_.getArray();
     119  Vector<Float> wei0 = weightSumNoMask_.getArray();
     120  Vector<uInt>  n = n_.getArray();
     121  Vector<uInt>  n0 = nNoMask_.getArray();
     122  cout << "S__" << "[" << spe[0] << "][" << spe[1] << "][" << spe[2] << "][" << spe[3] << "][" << spe[4] << "][" << spe[5] << "][" << spe[6] << "][" << spe[7] << "][" << spe[8] << "][" << spe[9] << "][" << spe[10] << "][" << spe[11] << "][" << spe[12] << "]" << endl;
     123  cout << "S0_" << "[" << spe0[0] << "][" << spe0[1] << "][" << spe0[2] << "][" << spe0[3] << "][" << spe0[4] << "][" << spe0[5] << "][" << spe0[6] << "][" << spe0[7] << "][" << spe0[8] << "][" << spe0[9] << "][" << spe0[10] << "][" << spe0[11] << "][" << spe0[12] << "]" << endl;
     124  cout << "W__" << "[" << wei[0] << "][" << wei[1] << "][" << wei[2] << "][" << wei[3] << "][" << wei[4] << "][" << wei[5] << "][" << wei[6] << "][" << wei[7] << "][" << wei[8] << "][" << wei[9] << "][" << wei[10] << "][" << wei[11] << "][" << wei[12] << "]" << endl;
     125  cout << "W0_" << "[" << wei0[0] << "][" << wei0[1] << "][" << wei0[2] << "][" << wei0[3] << "][" << wei0[4] << "][" << wei0[5] << "][" << wei0[6] << "][" << wei0[7] << "][" << wei0[8] << "][" << wei0[9] << "][" << wei0[10] << "][" << wei0[11] << "][" << wei0[12] << "]" << endl;
     126  cout << "N__" << "[" << n[0] << "][" << n[1] << "][" << n[2] << "][" << n[3] << "][" << n[4] << "][" << n[5] << "][" << n[6] << "][" << n[7] << "][" << n[8] << "][" << n[9] << "][" << n[10] << "][" << n[11] << "][" << n[12] << "]" << endl;
     127  cout << "N0_" << "[" << n0[0] << "][" << n0[1] << "][" << n0[2] << "][" << n0[3] << "][" << n0[4] << "][" << n0[5] << "][" << n0[6] << "][" << n0[7] << "][" << n0[8] << "][" << n0[9] << "][" << n0[10] << "][" << n0[11] << "][" << n0[12] << "]" << endl;
     128  cout << "***" << endl;
     129  //
     130}
     131
     132Float RowAccumulator::getTotalWeight(const MaskedArray<Float>& data,
     133                                     const Vector<Float>& tsys,
     134                                     const Double interval,
     135                                     const Double time,
     136                                     const Bool ignoreMask)
     137{
     138  Float totalWeight = 1.0;
     139
     140  Vector<Bool> m = data.getMask();
     141  if (!allEQ(m, False)) {  // only add these if not everything masked
     142    totalWeight *= addTsys(tsys, ignoreMask);
     143    totalWeight *= addInterval(interval, ignoreMask);
     144    addTime(time, ignoreMask);
     145  }
     146
     147  if (weightType_ == W_VAR) {
     148    Float fac = 1.0/variance(data);
     149    if (!ignoreMask && (m.nelements() == userMask_.nelements()))
     150      fac = 1.0/variance(data(userMask_));
     151
     152    totalWeight *= fac;
     153  }
     154
     155  return totalWeight;
     156}
     157
     158Float RowAccumulator::addTsys(const Vector<Float>& v, Bool ignoreMask)
    87159{
    88160  // @fixme this assume tsys is the same for all channels
    89161
    90162  Float w = 1.0;
    91   tsysSum_ += v[0];
    92   if ( weightType_ == asap::W_TSYS  || weightType_ == asap::W_TINTSYS ) {
     163  if (ignoreMask) {
     164    tsysSumNoMask_ += v[0];
     165  } else {
     166    tsysSum_ += v[0];
     167  }
     168  if ( weightType_ == W_TSYS  || weightType_ == W_TINTSYS ) {
    93169    w /= (v[0]*v[0]);
    94170  }
    … …  
    96172}
    97173
    98 void asap::RowAccumulator::addTime( casa::Double t )
    99 {
    100   timeSum_ += t;
    101 }
    102 
    103 Float asap::RowAccumulator::addInterval( casa::Double inter )
     174void RowAccumulator::addTime(Double t, Bool ignoreMask)
     175{
     176  if (ignoreMask) {
     177    timeSumNoMask_ += t;
     178  } else {
     179    timeSum_ += t;
     180  }
     181}
     182
     183Float RowAccumulator::addInterval(Double inter, Bool ignoreMask)
    104184{
    105185  Float w = 1.0;
    106   intervalSum_ += inter;
    107   if ( weightType_ == asap::W_TINT || weightType_ == asap::W_TINTSYS ) {
     186  if (ignoreMask) {
     187    intervalSumNoMask_ += inter;
     188  } else {
     189    intervalSum_ += inter;
     190  }
     191  if (weightType_ == W_TINT || weightType_ == W_TINTSYS) {
    108192    w /= Float(inter);
    109193  }
    … …  
    111195}
    112196
    113 void asap::RowAccumulator::reset( )
    114 {
    115   initialized_ = False;
    116   intervalSum_ = 0.0;
    117   tsysSum_.resize();
    118   timeSum_ = 0.0;
    119 }
    120 
    121 casa::Vector< casa::Float > RowAccumulator::getSpectrum( ) const
     197Vector<Float> RowAccumulator::getSpectrum() const
    122198{
    123199  return (spectrum_/weightSum_).getArray();
    124200}
    125201
    126 casa::Double asap::RowAccumulator::getTime( ) const
    127 {
    128   uInt n = max(n_);
    129   return timeSum_/Float(n);
    130 }
    131 
    132 casa::Double asap::RowAccumulator::getInterval( ) const
     202Double RowAccumulator::getTime() const
     203{
     204  return timeSum_/Float(max(n_));
     205}
     206
     207Double RowAccumulator::getInterval() const
    133208{
    134209  return intervalSum_;
    135210}
    136211
    137 casa::Vector< casa::Bool > RowAccumulator::getMask( ) const
     212Vector<Bool> RowAccumulator::getMask() const
    138213{
    139214  // Return the "total" mask - False where no points have been accumulated.
    … …  
    141216}
    142217
    143 casa::Vector< casa::Float > asap::RowAccumulator::getTsys( ) const
     218Vector<Float> RowAccumulator::getTsys() const
    144219{
    145220  // @fixme this assumes tsys.nelements() == 1
    … …  
    147222}
    148223
    149 void asap::RowAccumulator::setUserMask( const casa::Vector< casa::Bool > & m )
     224void RowAccumulator::setUserMask(const Vector<Bool>& m)
    150225{
    151226  userMask_.resize();
    … …  
    154229
    155230// Added by TT  check the state of RowAccumulator
    156 casa::Bool RowAccumulator::state() const
     231Bool RowAccumulator::state() const
    157232{
    158233  return initialized_;
    159234}
    160235
     236void RowAccumulator::replaceNaN()
     237{
     238  Vector<Float> v = spectrum_.getArray();
     239  Vector<Float> w = weightSum_.getArray();
     240  Vector<Float> vRef = spectrumNoMask_.getArray();
     241  Vector<Float> wRef = weightSumNoMask_.getArray();
     242
     243  //-------
     244  cout << "SB-";
     245  for (uInt i=0; i<13; ++i) {
     246    cout << "[" << v[i] << "]";
     247  }
     248  cout << endl;
     249  cout << "WB-";
     250  for (uInt i=0; i<13; ++i) {
     251    cout << "[" << w[i] << "]";
     252  }
     253  cout << endl;
     254  //-------
     255
     256  for (uInt i = 0; i < v.size(); ++i) {
     257    if (w[i] == 0.0) {
     258      v[i] = vRef[i];
     259      w[i] = wRef[i];
     260    }
     261  }
     262
     263  spectrum_.setData(v, Vector<Bool>(v.nelements(), True));
     264  weightSum_.setData(w, Vector<Bool>(w.nelements(), True));
     265
     266  //-------
     267  cout << "SA-";
     268  for (uInt i=0; i<13; ++i) {
     269    cout << "[" << v[i] << "]";
     270  }
     271  cout << endl;
     272  cout << "WA-";
     273  for (uInt i=0; i<13; ++i) {
     274    cout << "[" << w[i] << "]";
     275  }
     276  cout << endl;
     277  //-------
     278}

  • trunk/src/RowAccumulator.h

    r1819 r2125  
    1313#define ASAPROWACCUMULATOR_H
    1414
     15#include <math.h>
    1516#include <casa/aips.h>
    1617#include <casa/Arrays/Vector.h>
    … …  
    4849           const casa::Vector<casa::Bool>& m,
    4950           const casa::Vector<casa::Float>& tsys,
    50            casa::Double interval,
    51            casa::Double time);
     51           const casa::Double interval,
     52           const casa::Double time);
    5253  /**
    5354    * Also set a user mask which get combined with the individual masks
    … …  
    8485    * Reset the acummulator to the state at construction.
    8586    */
    86   void reset();
     87  void reset(const casa::uInt size=0, const casa::uInt tsysSize=0);
     88  void initialize(const casa::uInt size, const casa::uInt tsysSize);
    8789  /**
    8890    * check the initialization state
    8991    */
    9092  casa::Bool state() const;
     93  /**
     94    * replace NaN values with (normal) values at the same channels in the given spetrum.
     95    * (CAS-2776; 2011/04/07 by Wataru Kawasaki)
     96    */
     97  void replaceNaN();
    9198
    9299private:
    93   void addSpectrum( const casa::Vector<casa::Float>& v,
    94                     const casa::Vector<casa::Bool>& m,
    95                     casa::Float weight);
    96 
    97   casa::Float addTsys(const casa::Vector<casa::Float>& v);
    98   casa::Float addInterval(casa::Double inter);
    99   void addTime(casa::Double t);
     100  void addSpectrum(const casa::Vector<casa::Float>& v,
     101                   const casa::Vector<casa::Bool>& m,
     102                   const casa::Vector<casa::Float>& tsys,
     103                   const casa::Double interval,
     104                   const casa::Double time);
     105  void doAddSpectrum(const casa::Vector<casa::Float>& v,
     106                     const casa::Vector<casa::Bool>& m,
     107                     const casa::Vector<casa::Float>& tsys,
     108                     const casa::Double interval,
     109                     const casa::Double time,
     110                     const casa::Bool ignoreMask);
     111  casa::Float getTotalWeight(const casa::MaskedArray<casa::Float>& data,
     112                             const casa::Vector<casa::Float>& tsys,
     113                             const casa::Double interval,
     114                             const casa::Double time,
     115                             const casa::Bool ignoreMask);
     116  casa::Float addTsys(const casa::Vector<casa::Float>& v, casa::Bool ignoreMask);
     117  casa::Float addInterval(casa::Double inter, casa::Bool ignoreMask);
     118  void addTime(casa::Double t, casa::Bool ignoreMask);
    100119
    101120  WeightType weightType_;
    … …  
    106125  casa::MaskedArray<casa::uInt> n_;
    107126
     127  //these three are used for normalise() (CAS-2776; 2011/04/07 by WK)
     128  casa::MaskedArray<casa::Float> spectrumNoMask_;
     129  casa::MaskedArray<casa::Float> weightSumNoMask_;
     130  casa::MaskedArray<casa::uInt> nNoMask_;
     131
    108132  casa::Vector<casa::Bool> userMask_;
    109133
    110   casa::Vector<casa::Float> tsysSum_;
    111   casa::Double timeSum_;
    112   casa::Double intervalSum_;
     134  casa::Vector<casa::Float> tsysSum_, tsysSumNoMask_;
     135  casa::Double timeSum_, timeSumNoMask_;
     136  casa::Double intervalSum_, intervalSumNoMask_;
    113137};
    114138

  • trunk/src/STMath.cpp

    r2121 r2125  
    207207      ROScalarColumn<Double> tmp(tin, "TIME");
    208208      Double td;tmp.get(0,td);
    209       Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
    210                        && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
    211                        && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
     209      Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
     210                         && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
     211                         && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
    212212      Table subt;
    213213      if ( avmode == "SOURCE") {
    214         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
     214        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
    215215      } else if (avmode == "SCAN") {
    216         //subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
    217         subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO"))
    218                          && basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
     216        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
     217                      && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
    219218      } else {
    220219        subt = basesubt;
    … …  
    266265        acc.add(spec, !bflag, tsys, inter, time);
    267266      }
     267
     268      // If there exists a channel at which all the input spectra are masked,
     269      // spec has 'nan' values for that channel and it may affect the following
     270      // processes. To avoid this, replacing 'nan' values in spec with
     271      // weighted-mean of all spectra in the following line.
     272      // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
     273      acc.replaceNaN();
    268274    }
    269275    const Vector<Bool>& msk = acc.getMask();
    … …  
    278284      Vector<uChar> flg(msk.shape());
    279285      convertArray(flg, !msk);
     286      for (uInt k = 0; k < flg.nelements(); ++k) {
     287        uChar userFlag = 1 << 7;
     288        if (msk[k]==True) userFlag = 0 << 7;
     289        flg(k) = userFlag;
     290      }
     291
    280292      flagColOut.put(i, flg);
    281293      specColOut.put(i, acc.getSpectrum());
    … …  
    627639  // check if vector size is equal to nrow
    628640  Vector<Float> fact( val ) ;
    629   if ( fact.nelements() != in->nrow() ) {
     641  if (fact.nelements() != uInt(in->nrow())) {
    630642    throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
    631643  }
    … …  
    688700  // some checks
    689701  vector<uInt> nchans;
    690   for ( uInt i = 0 ; i < in->nrow() ; i++ ) {
    691     nchans.push_back( (in->getSpectrum( i )).size() ) ;
     702  for (Int i = 0 ; i < in->nrow() ; i++) {
     703    nchans.push_back((in->getSpectrum(i)).size());
    692704  }
    693705  //Vector<uInt> mchans( nchans ) ;
    … …  
    22482260  Table outtab = out->table();
    22492261
    2250   const uInt ntau = uInt(tau.size());
     2262  const Int ntau = uInt(tau.size());
    22512263  std::vector<float>::const_iterator tauit = tau.begin();
    22522264  AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
    … …  
    22922304{
    22932305  CountedPtr< Scantable > out = getScantable(in, false);
    2294   Table& table = out->table();
    2295   ArrayColumn<Float> specCol(table, "SPECTRA");
    2296   ArrayColumn<uChar> flagCol(table, "FLAGTRA");
    2297   Vector<Float> spec;
    2298   Vector<uChar> flag;
    2299   for ( uInt i=0; i<table.nrow(); ++i) {
    2300     specCol.get(i, spec);
    2301     flagCol.get(i, flag);
    2302     Vector<Bool> mask(flag.nelements());
    2303     convertArray(mask, flag);
    2304     Vector<Float> specout;
    2305     Vector<Bool> maskout;
    2306     if ( kernel == "hanning" ) {
    2307       mathutil::hanning(specout, maskout, spec , !mask);
    2308       convertArray(flag, !maskout);
    2309     } else if (  kernel == "rmedian" ) {
    2310       mathutil::runningMedian(specout, maskout, spec , mask, width);
    2311       convertArray(flag, maskout);
    2312     } else if ( kernel == "poly" ) {
    2313       mathutil::polyfit(specout, maskout, spec, !mask, width, order);
    2314       convertArray(flag, !maskout);
    2315     }
    2316     flagCol.put(i, flag);
    2317     specCol.put(i, specout);
     2306  Table table = out->table();
     2307
     2308  TableIterator iter(table, "IFNO");
     2309  while (!iter.pastEnd()) {
     2310    Table tab = iter.table();
     2311    ArrayColumn<Float> specCol(tab, "SPECTRA");
     2312    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     2313    Vector<Float> spec;
     2314    Vector<uChar> flag;
     2315    for (uInt i = 0; i < tab.nrow(); ++i) {
     2316      specCol.get(i, spec);
     2317      flagCol.get(i, flag);
     2318      Vector<Bool> mask(flag.nelements());
     2319      convertArray(mask, flag);
     2320      Vector<Float> specout;
     2321      Vector<Bool> maskout;
     2322      if (kernel == "hanning") {
     2323        mathutil::hanning(specout, maskout, spec, !mask);
     2324        convertArray(flag, !maskout);
     2325      } else if (kernel == "rmedian") {
     2326        mathutil::runningMedian(specout, maskout, spec , mask, width);
     2327        convertArray(flag, maskout);
     2328      } else if (kernel == "poly") {
     2329        mathutil::polyfit(specout, maskout, spec, !mask, width, order);
     2330        convertArray(flag, !maskout);
     2331      }
     2332
     2333      for (uInt j = 0; j < flag.nelements(); ++j) {
     2334        uChar userFlag = 1 << 7;
     2335        if (maskout[j]==True) userFlag = 0 << 7;
     2336        flag(j) = userFlag;
     2337      }
     2338
     2339      flagCol.put(i, flag);
     2340      specCol.put(i, specout);
     2341    }
     2342  ++iter;
    23182343  }
    23192344  return out;
    … …  
    28222847        if (fstart-1 > 0 ) {
    28232848          interp = spec[fstart-1];
    2824           if (fend+1 < spec.nelements()) {
     2849          if (fend+1 < Int(spec.nelements())) {
    28252850            interp = (interp+spec[fend+1])/2.0;
    28262851          }

  • trunk/src/python_Scantable.cpp

    r2111 r2125  
    155155    .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
    156156         (boost::python::arg("whichrow")=0) )
    157     //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline,
    158     // (boost::python::arg("thres")=3.0,
    159     //  boost::python::arg("niter")=1) )
     157    //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline)
    160158    //.def("_test_cin", &ScantableWrapper::testCin)
    161159  ;
Note: See TracChangeset for help on using the changeset viewer.