Ignore:
Timestamp:
03/30/11 21:43:43 (14 years ago)
Author:
WataruKawasaki
Message:

merged bug fixes from trunk (r2081)

Location:
branches/casa-prerelease/pre-asap
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/casa-prerelease/pre-asap

    • Property svn:mergeinfo changed
      /trunkmerged: 2081

  • branches/casa-prerelease/pre-asap/src

  • branches/casa-prerelease/pre-asap/src/SConscript

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/src/STLineFinder.cpp

    r2012 r2082  
    12031203    vel = scan->getAbcissa(last_row_used);
    12041204  } else {
    1205     for (int i = 0; i < spectrum.nelements(); ++i)
     1205    for (uInt i = 0; i < spectrum.nelements(); ++i)
    12061206      vel.push_back((double)i);
    12071207  }

  • branches/casa-prerelease/pre-asap/src/Scantable.cpp

    r2065 r2082  
    17641764{
    17651765  ofstream ofs;
    1766   String coordInfo;
     1766  String coordInfo = "";
    17671767  bool hasSameNchan = true;
    17681768  bool outTextFile = false;
    … …  
    17811781  Fitter fitter = Fitter();
    17821782  fitter.setExpression("poly", order);
     1783  //fitter.setIterClipping(thresClip, nIterClip);
    17831784
    17841785  int nRow = nrow();
    … …  
    17981799{
    17991800  ofstream ofs;
    1800   String coordInfo;
     1801  String coordInfo = "";
    18011802  bool hasSameNchan = true;
    18021803  bool outTextFile = false;
    … …  
    18151816  Fitter fitter = Fitter();
    18161817  fitter.setExpression("poly", order);
     1818  //fitter.setIterClipping(thresClip, nIterClip);
    18171819
    18181820  int nRow = nrow();
    … …  
    18561858}
    18571859
    1858 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1860void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile)
     1861{
    18591862  ofstream ofs;
    1860   String coordInfo;
     1863  String coordInfo = "";
    18611864  bool hasSameNchan = true;
    18621865  bool outTextFile = false;
    … …  
    18751878  //Fitter fitter = Fitter();
    18761879  //fitter.setExpression("cspline", nPiece);
     1880  //fitter.setIterClipping(thresClip, nIterClip);
    18771881
    18781882  int nRow = nrow();
    … …  
    18811885  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    18821886    chanMask = getCompositeChanMask(whichrow, mask);
    1883     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1887    //fitBaseline(chanMask, whichrow, fitter);
    18841888    //setSpectrum(fitter.getResidual(), whichrow);
    1885     std::vector<int> pieceRanges;
     1889    std::vector<int> pieceEdges;
    18861890    std::vector<float> params;
    1887     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1891    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    18881892    setSpectrum(res, whichrow);
    18891893    //
    18901894
    1891     std::vector<bool>  fixed;
    1892     fixed.clear();
    1893     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
     1895    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
    18941896  }
    18951897
    … …  
    19001902{
    19011903  ofstream ofs;
    1902   String coordInfo;
     1904  String coordInfo = "";
    19031905  bool hasSameNchan = true;
    19041906  bool outTextFile = false;
    … …  
    19171919  //Fitter fitter = Fitter();
    19181920  //fitter.setExpression("cspline", nPiece);
     1921  //fitter.setIterClipping(thresClip, nIterClip);
    19191922
    19201923  int nRow = nrow();
    … …  
    19501953
    19511954
    1952     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1955    //fitBaseline(chanMask, whichrow, fitter);
    19531956    //setSpectrum(fitter.getResidual(), whichrow);
    1954     std::vector<int> pieceRanges;
     1957    std::vector<int> pieceEdges;
    19551958    std::vector<float> params;
    1956     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1959    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19571960    setSpectrum(res, whichrow);
    19581961    //
    19591962
    1960     std::vector<bool> fixed;
    1961     fixed.clear();
    1962     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1963    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
    19631964  }
    19641965
    … …  
    19661967}
    19671968
    1968 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& idxEdge, std::vector<float>& params, int nPiece, float thresClip, int nIterClip, bool getResidual) {
     1969std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
     1970{
     1971  if (data.size() != mask.size()) {
     1972    throw(AipsError("data and mask sizes are not identical"));
     1973  }
    19691974  if (nPiece < 1) {
    19701975    throw(AipsError("wrong number of the sections for fitting"));
    1971   }
    1972   if (data.size() != mask.size()) {
    1973     throw(AipsError("data and mask have different size"));
    19741976  }
    19751977
    … …  
    19841986  }
    19851987
    1986   int nData = x.size();
    1987   int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
     1988  int initNData = x.size();
     1989
     1990  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    19881991  std::vector<double> invEdge;
    1989 
    19901992  idxEdge.clear();
    19911993  idxEdge.push_back(x[0]);
    1992 
    19931994  for (int i = 1; i < nPiece; ++i) {
    19941995    int valX = x[nElement*i];
    … …  
    19961997    invEdge.push_back(1.0/(double)valX);
    19971998  }
    1998 
    19991999  idxEdge.push_back(x[x.size()-1]+1);
    20002000
    2001   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2001  int nData = initNData;
     2002  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2003
     2004  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
    20022005  for (int i = 0; i < nChan; ++i) {
    20032006    double di = (double)i;
    … …  
    20112014    x3z1.push_back(dD*di*di*di);
    20122015    r1.push_back(0.0);
    2013   }
    2014 
    2015   int currentNData = nData;
    2016   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2016    residual.push_back(0.0);
     2017  }
    20172018
    20182019  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    … …  
    21282129      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    21292130        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2131        residual[i] = z1[i] - r1[i];
    21302132      }
    21312133      params.push_back(a0);
    … …  
    21472149      break;
    21482150    } else {
    2149       std::vector<double> rz;
    21502151      double stdDev = 0.0;
    21512152      for (int i = 0; i < nChan; ++i) {
    2152         double val = abs(r1[i] - z1[i]);
    2153         rz.push_back(val);
    2154         stdDev += val*val*(double)maskArray[i];
     2153        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    21552154      }
    21562155      stdDev = sqrt(stdDev/(double)nData);
    … …  
    21592158      int newNData = 0;
    21602159      for (int i = 0; i < nChan; ++i) {
    2161         if (rz[i] >= thres) {
     2160        if (abs(residual[i]) >= thres) {
    21622161          maskArray[i] = 0;
    21632162        }
    … …  
    21662165        }
    21672166      }
    2168       if (newNData == currentNData) {
     2167      if (newNData == nData) {
    21692168        break; //no more flag to add. iteration stops.
    21702169      } else {
    2171         currentNData = newNData;
     2170        nData = newNData;
    21722171      }
    21732172    }
    … …  
    21772176  if (getResidual) {
    21782177    for (int i = 0; i < nChan; ++i) {
    2179       result.push_back((float)(z1[i] - r1[i]));
     2178      result.push_back((float)residual[i]);
    21802179    }
    21812180  } else {
    … …  
    21882187}
    21892188
    2190 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     2189void 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)
     2190{
    21912191  ofstream ofs;
    2192   String coordInfo;
     2192  String coordInfo = "";
    21932193  bool hasSameNchan = true;
    21942194  bool outTextFile = false;
    … …  
    22062206
    22072207  //Fitter fitter = Fitter();
    2208   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2208  //fitter.setExpression("sinusoid", nWaves);
     2209  //fitter.setIterClipping(thresClip, nIterClip);
    22092210
    22102211  int nRow = nrow();
    … …  
    22132214  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22142215    chanMask = getCompositeChanMask(whichrow, mask);
    2215     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2216    //fitBaseline(chanMask, whichrow, fitter);
    22162217    //setSpectrum(fitter.getResidual(), whichrow);
    2217     std::vector<int> pieceRanges;
    22182218    std::vector<float> params;
    2219     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2219    //std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, true);
     2220    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22202221    setSpectrum(res, whichrow);
    22212222    //
    22222223
    2223     std::vector<bool>  fixed;
    2224     fixed.clear();
    2225     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
     2224    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
    22262225  }
    22272226
    … …  
    22292228}
    22302229
    2231 void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     2230void 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)
    22322231{
    22332232  ofstream ofs;
    2234   String coordInfo;
     2233  String coordInfo = "";
    22352234  bool hasSameNchan = true;
    22362235  bool outTextFile = false;
    … …  
    22482247
    22492248  //Fitter fitter = Fitter();
    2250   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2249  //fitter.setExpression("sinusoid", nWaves);
     2250  //fitter.setIterClipping(thresClip, nIterClip);
    22512251
    22522252  int nRow = nrow();
    … …  
    22822282
    22832283
    2284     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2284    //fitBaseline(chanMask, whichrow, fitter);
    22852285    //setSpectrum(fitter.getResidual(), whichrow);
    2286     std::vector<int> pieceRanges;
    22872286    std::vector<float> params;
    2288     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2287    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22892288    setSpectrum(res, whichrow);
    22902289    //
    22912290
    2292     std::vector<bool> fixed;
    2293     fixed.clear();
    2294     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
     2291    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
    22952292  }
    22962293
    … …  
    22982295}
    22992296
    2300 std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual) {
     2297std::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)
     2298{
    23012299  if (data.size() != mask.size()) {
    2302     throw(AipsError("data and mask have different size"));
    2303   }
    2304   if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
     2300    throw(AipsError("data and mask sizes are not identical"));
     2301  }
     2302  if (data.size() < 2) {
     2303    throw(AipsError("data size is too short"));
     2304  }
     2305  if (waveNumbers.size() == 0) {
     2306    throw(AipsError("missing wave number info"));
     2307  }
     2308  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
     2309  nWaves.reserve(waveNumbers.size());
     2310  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
     2311  sort(nWaves.begin(), nWaves.end());
     2312  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
     2313  nWaves.erase(end_it, nWaves.end());
     2314
     2315  int minNWaves = nWaves[0];
     2316  if (minNWaves < 0) {
    23052317    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
    2306   } else {
    2307     if (nMaxWavesInSW < nMinWavesInSW) {
    2308       int temp = nMaxWavesInSW;
    2309       nMaxWavesInSW = nMinWavesInSW;
    2310       nMinWavesInSW = temp;
    2311     }
    2312   }
     2318  }
     2319  bool hasConstantTerm = (minNWaves == 0);
    23132320
    23142321  int nChan = data.size();
    … …  
    23222329  }
    23232330
    2324   int nData = x.size();
    2325 
    2326   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2331  int initNData = x.size();
     2332
     2333  int nData = initNData;
     2334  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     2335
     2336  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
     2337  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.
     2338
     2339  // xArray : contains elemental values for computing the least-square matrix.
     2340  //          xArray.size() is nDOF and xArray[*].size() is nChan.
     2341  //          Each xArray element are as follows:
     2342  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
     2343  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
     2344  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
     2345  //          where (1 <= n <= nMaxWavesInSW),
     2346  //          or,
     2347  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
     2348  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
     2349  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
     2350  std::vector<std::vector<double> > xArray;
     2351  if (hasConstantTerm) {
     2352    std::vector<double> xu;
     2353    for (int j = 0; j < nChan; ++j) {
     2354      xu.push_back(1.0);
     2355    }
     2356    xArray.push_back(xu);
     2357  }
     2358  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
     2359    double xFactor = baseXFactor*(double)nWaves[i];
     2360    std::vector<double> xs, xc;
     2361    xs.clear();
     2362    xc.clear();
     2363    for (int j = 0; j < nChan; ++j) {
     2364      xs.push_back(sin(xFactor*(double)j));
     2365      xc.push_back(cos(xFactor*(double)j));
     2366    }
     2367    xArray.push_back(xs);
     2368    xArray.push_back(xc);
     2369  }
     2370
     2371  std::vector<double> z1, r1, residual;
    23272372  for (int i = 0; i < nChan; ++i) {
    2328     double di = (double)i;
    2329     double dD = (double)data[i];
    2330     x1.push_back(di);
    2331     x2.push_back(di*di);
    2332     x3.push_back(di*di*di);
    2333     z1.push_back(dD);
    2334     x1z1.push_back(di*dD);
    2335     x2z1.push_back(di*di*dD);
    2336     x3z1.push_back(di*di*di*dD);
     2373    z1.push_back((double)data[i]);
    23372374    r1.push_back(0.0);
    2338   }
    2339 
    2340   int currentNData = nData;
    2341   int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2375    residual.push_back(0.0);
     2376  }
    23422377
    23432378  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    2344     //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
    2345     //identity matrix (right).
    2346     //the right part is used to calculate the inverse matrix of the left part.
     2379    // xMatrix : horizontal concatenation of
     2380    //           the least-sq. matrix (left) and an
     2381    //           identity matrix (right).
     2382    // the right part is used to calculate the inverse matrix of the left part.
    23472383    double xMatrix[nDOF][2*nDOF];
    23482384    double zMatrix[nDOF];
    … …  
    23552391    }
    23562392
    2357     for (int i = 0; i < currentNData; ++i) {
    2358       if (maskArray[i] == 0) continue;
    2359       xMatrix[0][0] += 1.0;
    2360       xMatrix[0][1] += x1[i];
    2361       xMatrix[0][2] += x2[i];
    2362       xMatrix[0][3] += x3[i];
    2363       xMatrix[1][1] += x2[i];
    2364       xMatrix[1][2] += x3[i];
    2365       xMatrix[1][3] += x2[i]*x2[i];
    2366       xMatrix[2][2] += x2[i]*x2[i];
    2367       xMatrix[2][3] += x3[i]*x2[i];
    2368       xMatrix[3][3] += x3[i]*x3[i];
    2369       zMatrix[0] += z1[i];
    2370       zMatrix[1] += x1z1[i];
    2371       zMatrix[2] += x2z1[i];
    2372       zMatrix[3] += x3z1[i];
     2393    for (int k = 0; k < nChan; ++k) {
     2394      if (maskArray[k] == 0) continue;
     2395
     2396      for (int i = 0; i < nDOF; ++i) {
     2397        for (int j = i; j < nDOF; ++j) {
     2398          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
     2399        }
     2400        zMatrix[i] += z1[k] * xArray[i][k];
     2401      }
    23732402    }
    23742403
    … …  
    24112440    }
    24122441    //compute a vector y which consists of the coefficients of the sinusoids forming the
    2413     //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
    2414     //respectively.
     2442    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
     2443    //and cosine functions, respectively.
    24152444    std::vector<double> y;
     2445    params.clear();
    24162446    for (int i = 0; i < nDOF; ++i) {
    24172447      y.push_back(0.0);
    … …  
    24192449        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    24202450      }
    2421     }
    2422 
    2423     double a0 = y[0];
    2424     double a1 = y[1];
    2425     double a2 = y[2];
    2426     double a3 = y[3];
    2427     params.clear();
     2451      params.push_back(y[i]);
     2452    }
    24282453
    24292454    for (int i = 0; i < nChan; ++i) {
    2430       r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2431     }
    2432     params.push_back(a0);
    2433     params.push_back(a1);
    2434     params.push_back(a2);
    2435     params.push_back(a3);
    2436 
     2455      r1[i] = y[0];
     2456      for (int j = 1; j < nDOF; ++j) {
     2457        r1[i] += y[j]*xArray[j][i];
     2458      }
     2459      residual[i] = z1[i] - r1[i];
     2460    }
    24372461
    24382462    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    24392463      break;
    24402464    } else {
    2441       std::vector<double> rz;
    24422465      double stdDev = 0.0;
    24432466      for (int i = 0; i < nChan; ++i) {
    2444         double val = abs(r1[i] - z1[i]);
    2445         rz.push_back(val);
    2446         stdDev += val*val*(double)maskArray[i];
     2467        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    24472468      }
    24482469      stdDev = sqrt(stdDev/(double)nData);
    … …  
    24512472      int newNData = 0;
    24522473      for (int i = 0; i < nChan; ++i) {
    2453         if (rz[i] >= thres) {
     2474        if (abs(residual[i]) >= thres) {
    24542475          maskArray[i] = 0;
    24552476        }
    … …  
    24582479        }
    24592480      }
    2460       if (newNData == currentNData) {
    2461         break;                   //no additional flag. finish iteration.
     2481      if (newNData == nData) {
     2482        break; //no more flag to add. iteration stops.
    24622483      } else {
    2463         currentNData = newNData;
     2484        nData = newNData;
    24642485      }
    24652486    }
    … …  
    24692490  if (getResidual) {
    24702491    for (int i = 0; i < nChan; ++i) {
    2471       result.push_back((float)(z1[i] - r1[i]));
     2492      result.push_back((float)residual[i]);
    24722493    }
    24732494  } else {
    … …  
    24822503void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    24832504{
    2484   std::vector<double> abcsd = getAbcissa(whichrow);
    2485   std::vector<float> abcs;
    2486   for (uInt i = 0; i < abcsd.size(); ++i) {
    2487     abcs.push_back((float)abcsd[i]);
     2505  std::vector<double> dAbcissa = getAbcissa(whichrow);
     2506  std::vector<float> abcissa;
     2507  for (uInt i = 0; i < dAbcissa.size(); ++i) {
     2508    abcissa.push_back((float)dAbcissa[i]);
    24882509  }
    24892510  std::vector<float> spec = getSpectrum(whichrow);
    24902511
    2491   fitter.setData(abcs, spec, mask);
     2512  fitter.setData(abcissa, spec, mask);
    24922513  fitter.lfit();
    24932514}
    … …  
    25522573
    25532574/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2554 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, const std::vector<bool>& fixed) {
     2575void 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) {
    25552576  if (outLogger || outTextFile) {
    25562577    float rms = getRms(chanMask, whichrow);
    25572578    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2579    std::vector<bool> fixed;
     2580    fixed.clear();
    25582581
    25592582    if (outLogger) {
    … …  
    25682591
    25692592/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2570 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, const std::vector<bool>& fixed) {
     2593void 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) {
    25712594  if (outLogger || outTextFile) {
    25722595    float rms = getRms(chanMask, whichrow);
    25732596    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2597    std::vector<bool> fixed;
     2598    fixed.clear();
    25742599
    25752600    if (outLogger) {

  • branches/casa-prerelease/pre-asap/src/Scantable.h

    r2065 r2082  
    522522                               const std::string& blfile="");
    523523  void sinusoidBaseline(const std::vector<bool>& mask,
    524                         int nMinWavesInSW,
    525                         int nMaxWavesInSW,
     524                        const std::vector<int>& nWaves,
     525                        float maxWaveLength,
    526526                        float thresClip,
    527527                        int nIterClip,
     528                        bool getResidual=true,
    528529                        bool outLogger=false,
    529530                        const std::string& blfile="");
    530531  void autoSinusoidBaseline(const std::vector<bool>& mask,
    531                             int nMinWavesInSW,
    532                             int nMaxWavesInSW,
     532                            const std::vector<int>& nWaves,
     533                            float maxWaveLength,
    533534                            float thresClip,
    534535                            int nIterClip,
    … …  
    536537                            float threshold=3.0,
    537538                            int chanAvgLimit=1,
     539                            bool getResidual=true,
    538540                            bool outLogger=false,
    539541                            const std::string& blfile="");
    … …  
    672674  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    673675                                          const std::vector<bool>& mask,
     676                                          int nPiece,
    674677                                          std::vector<int>& idxEdge,
    675678                                          std::vector<float>& params,
    676                                           int nPiece=2,
    677679                                          float thresClip=3.0,
    678680                                          int nIterClip=1,
    … …  
    680682  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
    681683                                       const std::vector<bool>& mask,
    682                                        int nMinWavesInSW,
    683                                        int nMaxWavesInSW,
     684                                       const std::vector<int>& waveNumbers,
     685                                       float maxWaveLength,
    684686                                       std::vector<float>& params,
    685687                                       float thresClip=3.0,
    … …  
    698700  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    699701  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, Fitter& fitter);
    700   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, const std::vector<bool>& fixed);
    701   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, const std::vector<bool>& fixed);
     702  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);
     703  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);
    702704
    703705  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • branches/casa-prerelease/pre-asap/src/ScantableWrapper.h

    r2047 r2082  
    269269  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
    270270
    271   void sinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
    272   { table_->sinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile); }
    273 
    274   void autoSinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
    275   { table_->autoSinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
     271  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="")
     272  { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
     273
     274  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="")
     275  { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    276276
    277277  float getRms(const std::vector<bool>& mask, int whichrow)
Note: See TracChangeset for help on using the changeset viewer.