Changeset 2186 for trunk/src


Ignore:
Timestamp:
06/07/11 23:49:53 (13 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*sinusoid_baseline() params

Test Programs:

Put in Release Notes: Yes

Module(s):

Description: (1) Implemented an automated sinusoidal fitting functionality

(2) FFT available with scantable.fft()
(3) fixed a bug of parsing 'edge' param used by linefinder.
(4) a function to show progress status for row-based iterations.


Location:
trunk/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/MathUtils.cpp

    r2163 r2186  
    248248  }
    249249}
     250
     251void mathutil::doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
     252                                        std::vector<bool>& mask) {
     253  int fstart = -1;
     254  int fend = -1;
     255  for (uInt i = 0; i < mask.size(); ++i) {
     256    if (!mask[i]) {
     257      fstart = i;
     258      while (!mask[i] && i < mask.size()) {
     259        fend = i;
     260        i++;
     261      }
     262    }
     263
     264    // execute interpolation as the following criteria:
     265    // (1) for a masked region inside the spectrum, replace the spectral
     266    //     values with the mean of those at the two channels just outside
     267    //     the both edges of the masked region.
     268    // (2) for a masked region at the spectral edge, replace the values
     269    //     with the one at the nearest non-masked channel.
     270    //     (ZOH, but bilateral)
     271    Float interp = 0.0;
     272    if (fstart-1 > 0) {
     273      interp = data[fstart-1];
     274      if (fend+1 < Int(data.nelements())) {
     275        interp = (interp + data[fend+1]) / 2.0;
     276      }
     277    } else {
     278      interp = data[fend+1];
     279    }
     280    if (fstart > -1 && fend > -1) {
     281      for (int j = fstart; j <= fend; ++j) {
     282        data[j] = interp;
     283      }
     284    }
     285
     286    fstart = -1;
     287    fend = -1;
     288  }
     289}

  • trunk/src/MathUtils.h

    r1819 r2186  
    6868 * @param hwidth half-width of the smoothing window
    6969 */
    70  void runningMedian(casa::Vector<casa::Float>& out,
     70void runningMedian(casa::Vector<casa::Float>& out,
    7171                   casa::Vector<casa::Bool>& outflag,
    7272                   const casa::Vector<casa::Float>& in,
    … …  
    7474                   float hwidth);
    7575
    76  void polyfit(casa::Vector<casa::Float>& out,
    77                    casa::Vector<casa::Bool>& outmask,
    78                    const casa::Vector<casa::Float>& in,
    79                    const casa::Vector<casa::Bool>& mask,
    80                    float hwidth, int order);
     76void polyfit(casa::Vector<casa::Float>& out,
     77             casa::Vector<casa::Bool>& outmask,
     78             const casa::Vector<casa::Float>& in,
     79             const casa::Vector<casa::Bool>& mask,
     80             float hwidth, int order);
    8181
    8282// Generate specified statistic
    … …  
    8585
    8686// Return a position of min or max value
    87  casa::IPosition minMaxPos(const casa::String& which,
     87casa::IPosition minMaxPos(const casa::String& which,
    8888                 const casa::MaskedArray<casa::Float>& data);
    8989
    … …  
    106106casa::Vector<casa::String> toVectorString(const std::vector<std::string>& in);
    107107
     108void doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
     109                              std::vector<bool>& mask);
     110
    108111}
    109112

  • trunk/src/STMath.cpp

    r2177 r2186  
    4949#include <atnf/PKSIO/SrcType.h>
    5050
    51 #include "MathUtils.h"
    5251#include "RowAccumulator.h"
    5352#include "STAttr.h"
    … …  
    28272826  std::vector<float> res;
    28282827  Table tab = in->table();
    2829   ArrayColumn<Float> specCol(tab, "SPECTRA");
    2830   ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
    2831   FFTServer<Float,Complex> ffts;
     2828  std::vector<bool> mask;
    28322829
    28332830  if (whichrow.size() < 1) {          // for all rows (by default)
    28342831    int nrow = int(tab.nrow());
    28352832    for (int i = 0; i < nrow; ++i) {
    2836       Vector<Float> spec = specCol(i);
    2837       Vector<uChar> flag = flagCol(i);
    2838       doFFT(res, ffts, getRealImag, spec, flag);
     2833      res = in->execFFT(i, mask, getRealImag);
    28392834    }
    28402835  } else {                           // for specified rows
    28412836    for (uInt i = 0; i < whichrow.size(); ++i) {
    2842       Vector<Float> spec = specCol(whichrow[i]);
    2843       Vector<uChar> flag = flagCol(whichrow[i]);
    2844       doFFT(res, ffts, getRealImag, spec, flag);
     2837      res = in->execFFT(i, mask, getRealImag);
    28452838    }
    28462839  }
    … …  
    28492842}
    28502843
    2851 void asap::STMath::doFFT( std::vector<float>& out,
    2852                           FFTServer<Float,Complex>& ffts,
    2853                           bool getRealImag,
    2854                           Vector<Float>& spec,
    2855                           Vector<uChar>& flag )
    2856 {
    2857   doZeroOrderInterpolation(spec, flag);
    2858 
    2859   Vector<Complex> fftout;
    2860   ffts.fft0(fftout, spec);
    2861 
    2862   float norm = float(2.0/double(spec.size()));
    2863   if (getRealImag) {
    2864     for (uInt j = 0; j < fftout.size(); ++j) {
    2865       out.push_back(real(fftout[j])*norm);
    2866       out.push_back(imag(fftout[j])*norm);
    2867     }
    2868   } else {
    2869     for (uInt j = 0; j < fftout.size(); ++j) {
    2870       out.push_back(abs(fftout[j])*norm);
    2871       out.push_back(arg(fftout[j]));
    2872     }
    2873   }
    2874 
    2875 }
    2876 
    2877 void asap::STMath::doZeroOrderInterpolation( Vector<Float>& spec,
    2878                                              Vector<uChar>& flag )
    2879 {
    2880   int fstart = -1;
    2881   int fend = -1;
    2882   for (uInt i = 0; i < flag.nelements(); ++i) {
    2883     if (flag[i] > 0) {
    2884       fstart = i;
    2885       while (flag[i] > 0 && i < flag.nelements()) {
    2886         fend = i;
    2887         i++;
    2888       }
    2889     }
    2890 
    2891     // execute interpolation as the following criteria:
    2892     // (1) for a masked region inside the spectrum, replace the spectral
    2893     //     values with the mean of those at the two channels just outside
    2894     //     the both edges of the masked region.
    2895     // (2) for a masked region at the spectral edge, replace the values
    2896     //     with the one at the nearest non-masked channel.
    2897     //     (ZOH, but bilateral)
    2898     Float interp = 0.0;
    2899     if (fstart-1 > 0) {
    2900       interp = spec[fstart-1];
    2901       if (fend+1 < Int(spec.nelements())) {
    2902         interp = (interp + spec[fend+1]) / 2.0;
    2903       }
    2904     } else {
    2905       interp = spec[fend+1];
    2906     }
    2907     if (fstart > -1 && fend > -1) {
    2908       for (int j = fstart; j <= fend; ++j) {
    2909         spec[j] = interp;
    2910       }
    2911     }
    2912 
    2913     fstart = -1;
    2914     fend = -1;
    2915   }
    2916 }
    29172844
    29182845CountedPtr<Scantable>
    … …  
    29392866      Vector<Float> spec = specCol(i);
    29402867      Vector<uChar> flag = flagCol(i);
    2941       doZeroOrderInterpolation(spec, flag);
     2868      std::vector<bool> mask;
     2869      for (uInt j = 0; j < flag.nelements(); ++j) {
     2870        mask.push_back(!(flag[j]>0));
     2871      }
     2872      mathutil::doZeroOrderInterpolation(spec, mask);
    29422873
    29432874      Vector<Complex> lags;

  • trunk/src/STMath.h

    r2177 r2186  
    2121#include <casa/Utilities/CountedPtr.h>
    2222
    23 #include <scimath/Mathematics/FFTServer.h>
    2423#include <scimath/Mathematics/InterpolateArray1D.h>
    2524
    … …  
    433432                                        int index ) ;
    434433  double getMJD( string strtime ) ;
    435   void doFFT( std::vector<float>& out,
    436               casa::FFTServer< casa::Float, casa::Complex >& ffts,
    437               bool getRealImag,
    438               casa::Vector<casa::Float>& spec,
    439               casa::Vector<casa::uChar>& flag ) ;
    440   void doZeroOrderInterpolation( casa::Vector<casa::Float>& spec,
    441                                  casa::Vector<casa::uChar>& flag) ;
    442434
    443435  bool insitu_;

  • trunk/src/Scantable.cpp

    r2163 r2186  
    1111//
    1212#include <map>
    13 #include <fstream>
     13
     14#include <atnf/PKSIO/SrcType.h>
    1415
    1516#include <casa/aips.h>
     17#include <casa/iomanip.h>
    1618#include <casa/iostream.h>
    17 #include <casa/iomanip.h>
     19#include <casa/OS/File.h>
    1820#include <casa/OS/Path.h>
    19 #include <casa/OS/File.h>
    2021#include <casa/Arrays/Array.h>
     22#include <casa/Arrays/ArrayAccessor.h>
     23#include <casa/Arrays/ArrayLogical.h>
    2124#include <casa/Arrays/ArrayMath.h>
    2225#include <casa/Arrays/MaskArrMath.h>
    23 #include <casa/Arrays/ArrayLogical.h>
    24 #include <casa/Arrays/ArrayAccessor.h>
     26#include <casa/Arrays/Slice.h>
    2527#include <casa/Arrays/Vector.h>
    2628#include <casa/Arrays/VectorSTLIterator.h>
    27 #include <casa/Arrays/Slice.h>
    2829#include <casa/BasicMath/Math.h>
    2930#include <casa/BasicSL/Constants.h>
     31#include <casa/Containers/RecordField.h>
     32#include <casa/Logging/LogIO.h>
    3033#include <casa/Quanta/MVAngle.h>
    31 #include <casa/Containers/RecordField.h>
     34#include <casa/Quanta/MVTime.h>
    3235#include <casa/Utilities/GenSort.h>
    33 #include <casa/Logging/LogIO.h>
    34 
    35 #include <tables/Tables/TableParse.h>
    36 #include <tables/Tables/TableDesc.h>
    37 #include <tables/Tables/TableCopy.h>
    38 #include <tables/Tables/SetupNewTab.h>
    39 #include <tables/Tables/ScaColDesc.h>
    40 #include <tables/Tables/ArrColDesc.h>
    41 #include <tables/Tables/TableRow.h>
    42 #include <tables/Tables/TableVector.h>
    43 #include <tables/Tables/TableIter.h>
    44 
    45 #include <tables/Tables/ExprNode.h>
    46 #include <tables/Tables/TableRecord.h>
    47 #include <casa/Quanta/MVTime.h>
    48 #include <casa/Quanta/MVAngle.h>
    49 #include <measures/Measures/MeasRef.h>
    50 #include <measures/Measures/MeasTable.h>
     36
     37#include <coordinates/Coordinates/CoordinateUtil.h>
     38
    5139// needed to avoid error in .tcc
    5240#include <measures/Measures/MCDirection.h>
    5341//
    5442#include <measures/Measures/MDirection.h>
     43#include <measures/Measures/MEpoch.h>
    5544#include <measures/Measures/MFrequency.h>
    56 #include <measures/Measures/MEpoch.h>
     45#include <measures/Measures/MeasRef.h>
     46#include <measures/Measures/MeasTable.h>
     47#include <measures/TableMeasures/ScalarMeasColumn.h>
     48#include <measures/TableMeasures/TableMeasDesc.h>
    5749#include <measures/TableMeasures/TableMeasRefDesc.h>
    5850#include <measures/TableMeasures/TableMeasValueDesc.h>
    59 #include <measures/TableMeasures/TableMeasDesc.h>
    60 #include <measures/TableMeasures/ScalarMeasColumn.h>
    61 #include <coordinates/Coordinates/CoordinateUtil.h>
    62 
    63 #include <atnf/PKSIO/SrcType.h>
    64 #include "Scantable.h"
    65 #include "STPolLinear.h"
    66 #include "STPolCircular.h"
    67 #include "STPolStokes.h"
     51
     52#include <tables/Tables/ArrColDesc.h>
     53#include <tables/Tables/ExprNode.h>
     54#include <tables/Tables/ScaColDesc.h>
     55#include <tables/Tables/SetupNewTab.h>
     56#include <tables/Tables/TableCopy.h>
     57#include <tables/Tables/TableDesc.h>
     58#include <tables/Tables/TableIter.h>
     59#include <tables/Tables/TableParse.h>
     60#include <tables/Tables/TableRecord.h>
     61#include <tables/Tables/TableRow.h>
     62#include <tables/Tables/TableVector.h>
     63
     64#include "MathUtils.h"
    6865#include "STAttr.h"
    6966#include "STLineFinder.h"
    70 #include "MathUtils.h"
     67#include "STPolCircular.h"
     68#include "STPolLinear.h"
     69#include "STPolStokes.h"
     70#include "Scantable.h"
    7171
    7272using namespace casa;
    … …  
    18471847    setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    18481848    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
     1849    showProgressOnTerminal(whichrow, nRow);
    18491850  }
    18501851
    … …  
    19091910
    19101911    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     1912    showProgressOnTerminal(whichrow, nRow);
    19111913  }
    19121914
    … …  
    19501952
    19511953    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
     1954    showProgressOnTerminal(whichrow, nRow);
    19521955  }
    19531956
    … …  
    20182021
    20192022    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
     2023    showProgressOnTerminal(whichrow, nRow);
    20202024  }
    20212025
    … …  
    22432247}
    22442248
    2245 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
     2249  void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2250{
     2251  nWaves.clear();
     2252
     2253  if (applyFFT) {
     2254    string fftThAttr;
     2255    float fftThSigma;
     2256    int fftThTop;
     2257    parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
     2258    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
     2259  }
     2260
     2261  addAuxWaveNumbers(addNWaves, rejectNWaves, nWaves);
     2262}
     2263
     2264void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
     2265{
     2266  uInt idxSigma = fftThresh.find("sigma");
     2267  uInt idxTop   = fftThresh.find("top");
     2268
     2269  if (idxSigma == fftThresh.size() - 5) {
     2270    std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
     2271    is >> fftThSigma;
     2272    fftThAttr = "sigma";
     2273  } else if (idxTop == 0) {
     2274    std::istringstream is(fftThresh.substr(3));
     2275    is >> fftThTop;
     2276    fftThAttr = "top";
     2277  } else {
     2278    bool isNumber = true;
     2279    for (uInt i = 0; i < fftThresh.size()-1; ++i) {
     2280      char ch = (fftThresh.substr(i, 1).c_str())[0];
     2281      if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
     2282        isNumber = false;
     2283        break;
     2284      }
     2285    }
     2286    if (isNumber) {
     2287      std::istringstream is(fftThresh);
     2288      is >> fftThSigma;
     2289      fftThAttr = "sigma";
     2290    } else {
     2291      throw(AipsError("fftthresh has a wrong value"));
     2292    }
     2293  }
     2294}
     2295
     2296void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
     2297{
     2298  std::vector<float> fspec;
     2299  if (fftMethod == "fft") {
     2300    fspec = execFFT(whichrow, chanMask, false, true);
     2301  //} else if (fftMethod == "lsp") {
     2302  //  fspec = lombScarglePeriodogram(whichrow);
     2303  }
     2304
     2305  if (fftThAttr == "sigma") {
     2306    float mean  = 0.0;
     2307    float mean2 = 0.0;
     2308    for (uInt i = 0; i < fspec.size(); ++i) {
     2309      mean  += fspec[i];
     2310      mean2 += fspec[i]*fspec[i];
     2311    }
     2312    mean  /= float(fspec.size());
     2313    mean2 /= float(fspec.size());
     2314    float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
     2315
     2316    for (uInt i = 0; i < fspec.size(); ++i) {
     2317      if (fspec[i] >= thres) {
     2318        nWaves.push_back(i);
     2319      }
     2320    }
     2321
     2322  } else if (fftThAttr == "top") {
     2323    for (int i = 0; i < fftThTop; ++i) {
     2324      float max = 0.0;
     2325      int maxIdx = 0;
     2326      for (uInt j = 0; j < fspec.size(); ++j) {
     2327        if (fspec[j] > max) {
     2328          max = fspec[j];
     2329          maxIdx = j;
     2330        }
     2331      }
     2332      nWaves.push_back(maxIdx);
     2333      fspec[maxIdx] = 0.0;
     2334    }
     2335
     2336  }
     2337
     2338  if (nWaves.size() > 1) {
     2339    sort(nWaves.begin(), nWaves.end());
     2340  }
     2341}
     2342
     2343void Scantable::addAuxWaveNumbers(const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2344{
     2345  for (uInt i = 0; i < addNWaves.size(); ++i) {
     2346    bool found = false;
     2347    for (uInt j = 0; j < nWaves.size(); ++j) {
     2348      if (nWaves[j] == addNWaves[i]) {
     2349        found = true;
     2350        break;
     2351      }
     2352    }
     2353    if (!found) nWaves.push_back(addNWaves[i]);
     2354  }
     2355
     2356  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
     2357    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
     2358      if (*j == rejectNWaves[i]) {
     2359        j = nWaves.erase(j);
     2360      } else {
     2361        ++j;
     2362      }
     2363    }
     2364  }
     2365
     2366  if (nWaves.size() > 1) {
     2367    sort(nWaves.begin(), nWaves.end());
     2368    unique(nWaves.begin(), nWaves.end());
     2369  }
     2370}
     2371
     2372void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
    22462373{
    22472374  ofstream ofs;
    … …  
    22672394  int nRow = nrow();
    22682395  std::vector<bool> chanMask;
     2396  std::vector<int> nWaves;
    22692397
    22702398  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22712399    chanMask = getCompositeChanMask(whichrow, mask);
     2400    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
     2401
     2402    //FOR DEBUGGING------------
     2403    if (whichrow < 0) {// == nRow -1) {
     2404      cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
     2405      if (applyFFT) {
     2406          cout << "[ ";
     2407          for (uInt j = 0; j < nWaves.size(); ++j) {
     2408            cout << nWaves[j] << ", ";
     2409          }
     2410          cout << " ]    " << endl;
     2411      }
     2412      cout << flush;
     2413    }
     2414    //-------------------------
     2415
    22722416    //fitBaseline(chanMask, whichrow, fitter);
    22732417    //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    22742418    std::vector<float> params;
    2275     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     2419    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    22762420    setSpectrum(res, whichrow);
    22772421    //
    22782422
    22792423    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
     2424    showProgressOnTerminal(whichrow, nRow);
    22802425  }
    22812426
    … …  
    22832428}
    22842429
    2285 void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
     2430void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
    22862431{
    22872432  ofstream ofs;
    … …  
    23072452  int nRow = nrow();
    23082453  std::vector<bool> chanMask;
     2454  std::vector<int> nWaves;
     2455
    23092456  int minEdgeSize = getIFNos().size()*2;
    23102457  STLineFinder lineFinder = STLineFinder();
    … …  
    23362483    //-------------------------------------------------------
    23372484
     2485    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    23382486
    23392487    //fitBaseline(chanMask, whichrow, fitter);
    23402488    //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    23412489    std::vector<float> params;
    2342     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     2490    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    23432491    setSpectrum(res, whichrow);
    23442492    //
    23452493
    23462494    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
     2495    showProgressOnTerminal(whichrow, nRow);
    23472496  }
    23482497
    … …  
    23502499}
    23512500
    2352 std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, float maxWaveLength, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
     2501std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
    23532502{
    23542503  if (data.size() != mask.size()) {
    … …  
    23592508  }
    23602509  if (waveNumbers.size() == 0) {
    2361     throw(AipsError("missing wave number info"));
     2510    throw(AipsError("no wave numbers given"));
    23622511  }
    23632512  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
    … …  
    23902539
    23912540  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
    2392   double baseXFactor = 2.0*PI/(double)maxWaveLength/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
     2541  double baseXFactor = 2.0*PI/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
    23932542
    23942543  // xArray : contains elemental values for computing the least-square matrix.
    … …  
    25712720std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
    25722721{
    2573   std::vector<bool> chanMask = getMask(whichrow);
    2574   uInt chanMaskSize = chanMask.size();
    2575   if (chanMaskSize != inMask.size()) {
    2576     throw(AipsError("different mask sizes"));
    2577   }
    2578   for (uInt i = 0; i < chanMaskSize; ++i) {
    2579     chanMask[i] = chanMask[i] && inMask[i];
    2580   }
    2581 
    2582   return chanMask;
     2722  std::vector<bool> mask = getMask(whichrow);
     2723  uInt maskSize = mask.size();
     2724  if (maskSize != inMask.size()) {
     2725    throw(AipsError("mask sizes are not the same."));
     2726  }
     2727  for (uInt i = 0; i < maskSize; ++i) {
     2728    mask[i] = mask[i] && inMask[i];
     2729  }
     2730
     2731  return mask;
    25832732}
    25842733
    … …  
    26102759
    26112760/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
    2612 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter) {
     2761void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
     2762{
    26132763  if (outLogger || outTextFile) {
    26142764    std::vector<float> params = fitter.getParameters();
    … …  
    26282778
    26292779/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2630 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params) {
     2780void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params)
     2781{
    26312782  if (outLogger || outTextFile) {
    26322783    float rms = getRms(chanMask, whichrow);
    … …  
    26462797
    26472798/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2648 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params) {
     2799void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params)
     2800{
    26492801  if (outLogger || outTextFile) {
    26502802    float rms = getRms(chanMask, whichrow);
    … …  
    26632815}
    26642816
    2665 float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
     2817void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold)
     2818{
     2819  if (nTotal >= nTotalThreshold) {
     2820    int nInterval = int(floor(double(nTotal)/100.0));
     2821    if (nInterval == 0) nInterval++;
     2822
     2823    if (nProcessed == 0) {
     2824      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
     2825      printf("[  0%%]");
     2826      printf("\x1b[39m\x1b[0m");             //default attributes
     2827      fflush(NULL);
     2828    } else if (nProcessed % nInterval == 0) {
     2829      printf("\r\x1b[1C");                   //go to the 2nd column
     2830      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
     2831      printf("%3d", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
     2832      printf("\x1b[39m\x1b[0m");             //default attributes
     2833      printf("\x1b[2C");                     //go to the end of line
     2834      fflush(NULL);
     2835    }
     2836    if (nProcessed == nTotal - 1) {
     2837      printf("\r\x1b[K");                    //clear
     2838      fflush(NULL);
     2839    }
     2840  }
     2841}
     2842
     2843std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
     2844{
     2845  std::vector<bool>  mask = getMask(whichrow);
     2846
     2847  if (inMask.size() > 0) {
     2848    uInt maskSize = mask.size();
     2849    if (maskSize != inMask.size()) {
     2850      throw(AipsError("mask sizes are not the same."));
     2851    }
     2852    for (uInt i = 0; i < maskSize; ++i) {
     2853      mask[i] = mask[i] && inMask[i];
     2854    }
     2855  }
     2856
     2857  Vector<Float> spec = getSpectrum(whichrow);
     2858  mathutil::doZeroOrderInterpolation(spec, mask);
     2859
     2860  FFTServer<Float,Complex> ffts;
     2861  Vector<Complex> fftres;
     2862  ffts.fft0(fftres, spec);
     2863
     2864  std::vector<float> res;
     2865  float norm = float(2.0/double(spec.size()));
     2866
     2867  if (getRealImag) {
     2868    for (uInt i = 0; i < fftres.size(); ++i) {
     2869      res.push_back(real(fftres[i])*norm);
     2870      res.push_back(imag(fftres[i])*norm);
     2871    }
     2872  } else {
     2873    for (uInt i = 0; i < fftres.size(); ++i) {
     2874      res.push_back(abs(fftres[i])*norm);
     2875      if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
     2876    }
     2877  }
     2878
     2879  return res;
     2880}
     2881
     2882
     2883float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
     2884{
    26662885  Vector<Float> spec;
    26672886  specCol_.get(whichrow, spec);
    … …  
    26852904
    26862905
    2687 std::string Scantable::formatBaselineParamsHeader(int whichrow,
    2688                                                   const std::string& masklist,
    2689                                                   bool verbose) const
     2906std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
    26902907{
    26912908  ostringstream oss;
    … …  
    27222939}
    27232940
    2724   std::string Scantable::formatBaselineParams(const std::vector<float>& params,
    2725                                               const std::vector<bool>& fixed,
    2726                                               float rms,
    2727                                               const std::string& masklist,
    2728                                               int whichrow,
    2729                                               bool verbose,
    2730                                               int start, int count,
    2731                                               bool resetparamid) const
     2941std::string Scantable::formatBaselineParams(const std::vector<float>& params,
     2942                                            const std::vector<bool>& fixed,
     2943                                            float rms,
     2944                                            const std::string& masklist,
     2945                                            int whichrow,
     2946                                            bool verbose,
     2947                                            int start, int count,
     2948                                            bool resetparamid) const
    27322949{
    27332950  int nParam = (int)(params.size());
    … …  
    29043121}
    29053122
    2906 /*
    2907 STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
    2908 {
    2909   Fitter fitter = Fitter();
    2910   fitter.setExpression("poly", order);
    2911 
    2912   std::vector<bool> fmask = getMask(rowno);
    2913   if (fmask.size() != mask.size()) {
    2914     throw(AipsError("different mask sizes"));
    2915   }
    2916   for (int i = 0; i < fmask.size(); ++i) {
    2917     fmask[i] = fmask[i] && mask[i];
    2918   }
    2919 
    2920   fitBaseline(fmask, rowno, fitter);
    2921   setSpectrum(fitter.getResidual(), rowno);
    2922   return fitter.getFitEntry();
    2923 }
    2924 */
    29253123
    29263124}

  • trunk/src/Scantable.h

    r2163 r2186  
    1414
    1515// STL
     16#include <fstream>
     17#include <iostream>
     18#include <sstream>
    1619#include <string>
    1720#include <vector>
    18 #include <iostream>
    19 #include <fstream>
    2021// AIPS++
    2122#include <casa/aips.h>
     23#include <casa/Arrays/MaskedArray.h>
     24#include <casa/Arrays/Vector.h>
     25#include <casa/BasicSL/String.h>
    2226#include <casa/Containers/Record.h>
    23 #include <casa/Arrays/MaskedArray.h>
    24 #include <casa/BasicSL/String.h>
     27#include <casa/Exceptions/Error.h>
     28#include <casa/Quanta/Quantum.h>
    2529#include <casa/Utilities/CountedPtr.h>
    2630
    27 #include <tables/Tables/Table.h>
     31#include <coordinates/Coordinates/SpectralCoordinate.h>
     32
     33#include <measures/TableMeasures/ScalarMeasColumn.h>
     34
     35#include <scimath/Mathematics/FFTServer.h>
     36
    2837#include <tables/Tables/ArrayColumn.h>
    2938#include <tables/Tables/ScalarColumn.h>
    30 
    31 #include <measures/TableMeasures/ScalarMeasColumn.h>
    32 
    33 #include <coordinates/Coordinates/SpectralCoordinate.h>
    34 
    35 #include <casa/Arrays/Vector.h>
    36 #include <casa/Quanta/Quantum.h>
    37 
    38 #include <casa/Exceptions/Error.h>
     39#include <tables/Tables/Table.h>
    3940
    4041#include "Logger.h"
    41 #include "STHeader.h"
    42 #include "STFrequencies.h"
    43 #include "STWeather.h"
    44 #include "STFocus.h"
    45 #include "STTcal.h"
    46 #include "STMolecules.h"
    47 #include "STSelector.h"
    48 #include "STHistory.h"
    49 #include "STPol.h"
     42#include "MathUtils.h"
    5043#include "STFit.h"
    5144#include "STFitEntry.h"
    5245#include "STFitter.h"
     46#include "STFocus.h"
     47#include "STFrequencies.h"
     48#include "STHeader.h"
     49#include "STHistory.h"
     50#include "STMolecules.h"
     51#include "STPol.h"
     52#include "STSelector.h"
     53#include "STTcal.h"
     54#include "STWeather.h"
    5355
    5456namespace asap {
    … …  
    528530                               const std::string& blfile="");
    529531  void sinusoidBaseline(const std::vector<bool>& mask,
    530                         const std::vector<int>& nWaves,
    531                         float maxWaveLength,
     532                        const bool applyFFT,
     533                        const std::string& fftMethod,
     534                        const std::string& fftThresh,
     535                        const std::vector<int>& addNWaves,
     536                        const std::vector<int>& rejectNWaves,
    532537                        float thresClip,
    533538                        int nIterClip,
    … …  
    536541                        const std::string& blfile="");
    537542  void autoSinusoidBaseline(const std::vector<bool>& mask,
    538                             const std::vector<int>& nWaves,
    539                             float maxWaveLength,
     543                            const bool applyFFT,
     544                            const std::string& fftMethod,
     545                            const std::string& fftThresh,
     546                            const std::vector<int>& addNWaves,
     547                            const std::vector<int>& rejectNWaves,
    540548                            float thresClip,
    541549                            int nIterClip,
    … …  
    546554                            bool outLogger=false,
    547555                            const std::string& blfile="");
     556  std::vector<float> execFFT(const int whichrow,
     557                             const std::vector<bool>& inMask,
     558                             bool getRealImag=false,
     559                             bool getAmplitudeOnly=false);
    548560  float getRms(const std::vector<bool>& mask, int whichrow);
    549561  std::string formatBaselineParams(const std::vector<float>& params,
    … …  
    689701                                       const std::vector<bool>& mask,
    690702                                       const std::vector<int>& waveNumbers,
    691                                        float maxWaveLength,
    692703                                       std::vector<float>& params,
    693704                                       float thresClip=3.0,
    694705                                       int nIterClip=1,
    695706                                       bool getResidual=true);
     707  void selectWaveNumbers(const int whichrow,
     708                         const std::vector<bool>& chanMask,
     709                         const bool applyFFT,
     710                         const std::string& fftMethod,
     711                         const std::string& fftThresh,
     712                         const std::vector<int>& addNWaves,
     713                         const std::vector<int>& rejectNWaves,
     714                         std::vector<int>& nWaves);
     715  void parseThresholdExpression(const std::string& fftThresh,
     716                                std::string& fftThAttr,
     717                                float& fftThSigma,
     718                                int& fftThTop);
     719  void doSelectWaveNumbers(const int whichrow,
     720                           const std::vector<bool>& chanMask,
     721                           const std::string& fftMethod,
     722                           const float fftThSigma,
     723                           const int fftThTop,
     724                           const std::string& fftThAttr,
     725                           std::vector<int>& nWaves);
     726  void addAuxWaveNumbers(const std::vector<int>& addNWaves,
     727                         const std::vector<int>& rejectNWaves,
     728                         std::vector<int>& nWaves);
    696729  bool hasSameNchanOverIFs();
    697730  std::string getMaskRangeList(const std::vector<bool>& mask,
    … …  
    708741  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params);
    709742  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<float>& params);
     743  void showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold=1000);
    710744
    711745  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • trunk/src/ScantableWrapper.h

    r2178 r2186  
    276276  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    277277
    278   void sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
    279   { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
    280 
    281   void autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
    282   { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
     278  void sinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
     279  { table_->sinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, getresidual, outlog, blfile); }
     280
     281  void autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
     282  { table_->autoSinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    283283
    284284  float getRms(const std::vector<bool>& mask, int whichrow)
    … …  
    293293  bool getFlagtraFast(int whichrow=0) const
    294294  { return table_->getFlagtraFast(casa::uInt(whichrow)); }
     295
     296  std::vector<float> execFFT(int whichrow, const std::vector<bool>& mask, bool getRealImag=false, bool getAmplitudeOnly=false)
     297  { return table_->execFFT(whichrow, mask, getRealImag, getAmplitudeOnly); }
    295298
    296299

  • trunk/src/python_Scantable.cpp

    r2178 r2186  
    154154    .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
    155155         (boost::python::arg("whichrow")=0) )
     156    .def("_fft", &ScantableWrapper::execFFT)
    156157    //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline)
    157158    //.def("_test_cin", &ScantableWrapper::testCin)
Note: See TracChangeset for help on using the changeset viewer.