Changeset 2081 for trunk/src


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

New Development: No

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): Scantable

Description: Scantable::sinusoidBaseline(), Scantable::autoSinusoidBaseline()


Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STLineFinder.cpp

    r2012 r2081  
    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  }

  • trunk/src/Scantable.cpp

    r2080 r2081  
    17781778{
    17791779  ofstream ofs;
    1780   String coordInfo;
     1780  String coordInfo = "";
    17811781  bool hasSameNchan = true;
    17821782  bool outTextFile = false;
    … …  
    17951795  Fitter fitter = Fitter();
    17961796  fitter.setExpression("poly", order);
     1797  //fitter.setIterClipping(thresClip, nIterClip);
    17971798
    17981799  int nRow = nrow();
    … …  
    18121813{
    18131814  ofstream ofs;
    1814   String coordInfo;
     1815  String coordInfo = "";
    18151816  bool hasSameNchan = true;
    18161817  bool outTextFile = false;
    … …  
    18291830  Fitter fitter = Fitter();
    18301831  fitter.setExpression("poly", order);
     1832  //fitter.setIterClipping(thresClip, nIterClip);
    18311833
    18321834  int nRow = nrow();
    … …  
    18701872}
    18711873
    1872 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1874void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile)
     1875{
    18731876  ofstream ofs;
    1874   String coordInfo;
     1877  String coordInfo = "";
    18751878  bool hasSameNchan = true;
    18761879  bool outTextFile = false;
    … …  
    18891892  //Fitter fitter = Fitter();
    18901893  //fitter.setExpression("cspline", nPiece);
     1894  //fitter.setIterClipping(thresClip, nIterClip);
    18911895
    18921896  int nRow = nrow();
    … …  
    18951899  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    18961900    chanMask = getCompositeChanMask(whichrow, mask);
    1897     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1901    //fitBaseline(chanMask, whichrow, fitter);
    18981902    //setSpectrum(fitter.getResidual(), whichrow);
    1899     std::vector<int> pieceRanges;
     1903    std::vector<int> pieceEdges;
    19001904    std::vector<float> params;
    1901     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1905    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19021906    setSpectrum(res, whichrow);
    19031907    //
    19041908
    1905     std::vector<bool>  fixed;
    1906     fixed.clear();
    1907     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
     1909    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
    19081910  }
    19091911
    … …  
    19141916{
    19151917  ofstream ofs;
    1916   String coordInfo;
     1918  String coordInfo = "";
    19171919  bool hasSameNchan = true;
    19181920  bool outTextFile = false;
    … …  
    19311933  //Fitter fitter = Fitter();
    19321934  //fitter.setExpression("cspline", nPiece);
     1935  //fitter.setIterClipping(thresClip, nIterClip);
    19331936
    19341937  int nRow = nrow();
    … …  
    19641967
    19651968
    1966     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1969    //fitBaseline(chanMask, whichrow, fitter);
    19671970    //setSpectrum(fitter.getResidual(), whichrow);
    1968     std::vector<int> pieceRanges;
     1971    std::vector<int> pieceEdges;
    19691972    std::vector<float> params;
    1970     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1973    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19711974    setSpectrum(res, whichrow);
    19721975    //
    19731976
    1974     std::vector<bool> fixed;
    1975     fixed.clear();
    1976     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1977    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
    19771978  }
    19781979
    … …  
    19801981}
    19811982
    1982 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) {
     1983std::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)
     1984{
     1985  if (data.size() != mask.size()) {
     1986    throw(AipsError("data and mask sizes are not identical"));
     1987  }
    19831988  if (nPiece < 1) {
    19841989    throw(AipsError("wrong number of the sections for fitting"));
    1985   }
    1986   if (data.size() != mask.size()) {
    1987     throw(AipsError("data and mask have different size"));
    19881990  }
    19891991
    … …  
    19982000  }
    19992001
    2000   int nData = x.size();
    2001   int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
     2002  int initNData = x.size();
     2003
     2004  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    20022005  std::vector<double> invEdge;
    2003 
    20042006  idxEdge.clear();
    20052007  idxEdge.push_back(x[0]);
    2006 
    20072008  for (int i = 1; i < nPiece; ++i) {
    20082009    int valX = x[nElement*i];
    … …  
    20102011    invEdge.push_back(1.0/(double)valX);
    20112012  }
    2012 
    20132013  idxEdge.push_back(x[x.size()-1]+1);
    20142014
    2015   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2015  int nData = initNData;
     2016  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2017
     2018  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
    20162019  for (int i = 0; i < nChan; ++i) {
    20172020    double di = (double)i;
    … …  
    20252028    x3z1.push_back(dD*di*di*di);
    20262029    r1.push_back(0.0);
    2027   }
    2028 
    2029   int currentNData = nData;
    2030   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2030    residual.push_back(0.0);
     2031  }
    20312032
    20322033  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    … …  
    21422143      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    21432144        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2145        residual[i] = z1[i] - r1[i];
    21442146      }
    21452147      params.push_back(a0);
    … …  
    21612163      break;
    21622164    } else {
    2163       std::vector<double> rz;
    21642165      double stdDev = 0.0;
    21652166      for (int i = 0; i < nChan; ++i) {
    2166         double val = abs(r1[i] - z1[i]);
    2167         rz.push_back(val);
    2168         stdDev += val*val*(double)maskArray[i];
     2167        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    21692168      }
    21702169      stdDev = sqrt(stdDev/(double)nData);
    … …  
    21732172      int newNData = 0;
    21742173      for (int i = 0; i < nChan; ++i) {
    2175         if (rz[i] >= thres) {
     2174        if (abs(residual[i]) >= thres) {
    21762175          maskArray[i] = 0;
    21772176        }
    … …  
    21802179        }
    21812180      }
    2182       if (newNData == currentNData) {
     2181      if (newNData == nData) {
    21832182        break; //no more flag to add. iteration stops.
    21842183      } else {
    2185         currentNData = newNData;
     2184        nData = newNData;
    21862185      }
    21872186    }
    … …  
    21912190  if (getResidual) {
    21922191    for (int i = 0; i < nChan; ++i) {
    2193       result.push_back((float)(z1[i] - r1[i]));
     2192      result.push_back((float)residual[i]);
    21942193    }
    21952194  } else {
    … …  
    22022201}
    22032202
    2204 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     2203void 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)
     2204{
    22052205  ofstream ofs;
    2206   String coordInfo;
     2206  String coordInfo = "";
    22072207  bool hasSameNchan = true;
    22082208  bool outTextFile = false;
    … …  
    22202220
    22212221  //Fitter fitter = Fitter();
    2222   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2222  //fitter.setExpression("sinusoid", nWaves);
     2223  //fitter.setIterClipping(thresClip, nIterClip);
    22232224
    22242225  int nRow = nrow();
    … …  
    22272228  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22282229    chanMask = getCompositeChanMask(whichrow, mask);
    2229     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2230    //fitBaseline(chanMask, whichrow, fitter);
    22302231    //setSpectrum(fitter.getResidual(), whichrow);
    2231     std::vector<int> pieceRanges;
    22322232    std::vector<float> params;
    2233     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2233    //std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, true);
     2234    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22342235    setSpectrum(res, whichrow);
    22352236    //
    22362237
    2237     std::vector<bool>  fixed;
    2238     fixed.clear();
    2239     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
     2238    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
    22402239  }
    22412240
    … …  
    22432242}
    22442243
    2245 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)
     2244void 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)
    22462245{
    22472246  ofstream ofs;
    2248   String coordInfo;
     2247  String coordInfo = "";
    22492248  bool hasSameNchan = true;
    22502249  bool outTextFile = false;
    … …  
    22622261
    22632262  //Fitter fitter = Fitter();
    2264   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2263  //fitter.setExpression("sinusoid", nWaves);
     2264  //fitter.setIterClipping(thresClip, nIterClip);
    22652265
    22662266  int nRow = nrow();
    … …  
    22962296
    22972297
    2298     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2298    //fitBaseline(chanMask, whichrow, fitter);
    22992299    //setSpectrum(fitter.getResidual(), whichrow);
    2300     std::vector<int> pieceRanges;
    23012300    std::vector<float> params;
    2302     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2301    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    23032302    setSpectrum(res, whichrow);
    23042303    //
    23052304
    2306     std::vector<bool> fixed;
    2307     fixed.clear();
    2308     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
     2305    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
    23092306  }
    23102307
    … …  
    23122309}
    23132310
    2314 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) {
     2311std::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)
     2312{
    23152313  if (data.size() != mask.size()) {
    2316     throw(AipsError("data and mask have different size"));
    2317   }
    2318   if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
     2314    throw(AipsError("data and mask sizes are not identical"));
     2315  }
     2316  if (data.size() < 2) {
     2317    throw(AipsError("data size is too short"));
     2318  }
     2319  if (waveNumbers.size() == 0) {
     2320    throw(AipsError("missing wave number info"));
     2321  }
     2322  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
     2323  nWaves.reserve(waveNumbers.size());
     2324  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
     2325  sort(nWaves.begin(), nWaves.end());
     2326  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
     2327  nWaves.erase(end_it, nWaves.end());
     2328
     2329  int minNWaves = nWaves[0];
     2330  if (minNWaves < 0) {
    23192331    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
    2320   } else {
    2321     if (nMaxWavesInSW < nMinWavesInSW) {
    2322       int temp = nMaxWavesInSW;
    2323       nMaxWavesInSW = nMinWavesInSW;
    2324       nMinWavesInSW = temp;
    2325     }
    2326   }
     2332  }
     2333  bool hasConstantTerm = (minNWaves == 0);
    23272334
    23282335  int nChan = data.size();
    … …  
    23362343  }
    23372344
    2338   int nData = x.size();
    2339 
    2340   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2345  int initNData = x.size();
     2346
     2347  int nData = initNData;
     2348  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     2349
     2350  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
     2351  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.
     2352
     2353  // xArray : contains elemental values for computing the least-square matrix.
     2354  //          xArray.size() is nDOF and xArray[*].size() is nChan.
     2355  //          Each xArray element are as follows:
     2356  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
     2357  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
     2358  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
     2359  //          where (1 <= n <= nMaxWavesInSW),
     2360  //          or,
     2361  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
     2362  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
     2363  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
     2364  std::vector<std::vector<double> > xArray;
     2365  if (hasConstantTerm) {
     2366    std::vector<double> xu;
     2367    for (int j = 0; j < nChan; ++j) {
     2368      xu.push_back(1.0);
     2369    }
     2370    xArray.push_back(xu);
     2371  }
     2372  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
     2373    double xFactor = baseXFactor*(double)nWaves[i];
     2374    std::vector<double> xs, xc;
     2375    xs.clear();
     2376    xc.clear();
     2377    for (int j = 0; j < nChan; ++j) {
     2378      xs.push_back(sin(xFactor*(double)j));
     2379      xc.push_back(cos(xFactor*(double)j));
     2380    }
     2381    xArray.push_back(xs);
     2382    xArray.push_back(xc);
     2383  }
     2384
     2385  std::vector<double> z1, r1, residual;
    23412386  for (int i = 0; i < nChan; ++i) {
    2342     double di = (double)i;
    2343     double dD = (double)data[i];
    2344     x1.push_back(di);
    2345     x2.push_back(di*di);
    2346     x3.push_back(di*di*di);
    2347     z1.push_back(dD);
    2348     x1z1.push_back(di*dD);
    2349     x2z1.push_back(di*di*dD);
    2350     x3z1.push_back(di*di*di*dD);
     2387    z1.push_back((double)data[i]);
    23512388    r1.push_back(0.0);
    2352   }
    2353 
    2354   int currentNData = nData;
    2355   int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2389    residual.push_back(0.0);
     2390  }
    23562391
    23572392  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    2358     //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
    2359     //identity matrix (right).
    2360     //the right part is used to calculate the inverse matrix of the left part.
     2393    // xMatrix : horizontal concatenation of
     2394    //           the least-sq. matrix (left) and an
     2395    //           identity matrix (right).
     2396    // the right part is used to calculate the inverse matrix of the left part.
    23612397    double xMatrix[nDOF][2*nDOF];
    23622398    double zMatrix[nDOF];
    … …  
    23692405    }
    23702406
    2371     for (int i = 0; i < currentNData; ++i) {
    2372       if (maskArray[i] == 0) continue;
    2373       xMatrix[0][0] += 1.0;
    2374       xMatrix[0][1] += x1[i];
    2375       xMatrix[0][2] += x2[i];
    2376       xMatrix[0][3] += x3[i];
    2377       xMatrix[1][1] += x2[i];
    2378       xMatrix[1][2] += x3[i];
    2379       xMatrix[1][3] += x2[i]*x2[i];
    2380       xMatrix[2][2] += x2[i]*x2[i];
    2381       xMatrix[2][3] += x3[i]*x2[i];
    2382       xMatrix[3][3] += x3[i]*x3[i];
    2383       zMatrix[0] += z1[i];
    2384       zMatrix[1] += x1z1[i];
    2385       zMatrix[2] += x2z1[i];
    2386       zMatrix[3] += x3z1[i];
     2407    for (int k = 0; k < nChan; ++k) {
     2408      if (maskArray[k] == 0) continue;
     2409
     2410      for (int i = 0; i < nDOF; ++i) {
     2411        for (int j = i; j < nDOF; ++j) {
     2412          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
     2413        }
     2414        zMatrix[i] += z1[k] * xArray[i][k];
     2415      }
    23872416    }
    23882417
    … …  
    24252454    }
    24262455    //compute a vector y which consists of the coefficients of the sinusoids forming the
    2427     //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
    2428     //respectively.
     2456    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
     2457    //and cosine functions, respectively.
    24292458    std::vector<double> y;
     2459    params.clear();
    24302460    for (int i = 0; i < nDOF; ++i) {
    24312461      y.push_back(0.0);
    … …  
    24332463        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    24342464      }
    2435     }
    2436 
    2437     double a0 = y[0];
    2438     double a1 = y[1];
    2439     double a2 = y[2];
    2440     double a3 = y[3];
    2441     params.clear();
     2465      params.push_back(y[i]);
     2466    }
    24422467
    24432468    for (int i = 0; i < nChan; ++i) {
    2444       r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2445     }
    2446     params.push_back(a0);
    2447     params.push_back(a1);
    2448     params.push_back(a2);
    2449     params.push_back(a3);
    2450 
     2469      r1[i] = y[0];
     2470      for (int j = 1; j < nDOF; ++j) {
     2471        r1[i] += y[j]*xArray[j][i];
     2472      }
     2473      residual[i] = z1[i] - r1[i];
     2474    }
    24512475
    24522476    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    24532477      break;
    24542478    } else {
    2455       std::vector<double> rz;
    24562479      double stdDev = 0.0;
    24572480      for (int i = 0; i < nChan; ++i) {
    2458         double val = abs(r1[i] - z1[i]);
    2459         rz.push_back(val);
    2460         stdDev += val*val*(double)maskArray[i];
     2481        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    24612482      }
    24622483      stdDev = sqrt(stdDev/(double)nData);
    … …  
    24652486      int newNData = 0;
    24662487      for (int i = 0; i < nChan; ++i) {
    2467         if (rz[i] >= thres) {
     2488        if (abs(residual[i]) >= thres) {
    24682489          maskArray[i] = 0;
    24692490        }
    … …  
    24722493        }
    24732494      }
    2474       if (newNData == currentNData) {
    2475         break;                   //no additional flag. finish iteration.
     2495      if (newNData == nData) {
     2496        break; //no more flag to add. iteration stops.
    24762497      } else {
    2477         currentNData = newNData;
     2498        nData = newNData;
    24782499      }
    24792500    }
    … …  
    24832504  if (getResidual) {
    24842505    for (int i = 0; i < nChan; ++i) {
    2485       result.push_back((float)(z1[i] - r1[i]));
     2506      result.push_back((float)residual[i]);
    24862507    }
    24872508  } else {
    … …  
    24962517void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    24972518{
    2498   std::vector<double> abcsd = getAbcissa(whichrow);
    2499   std::vector<float> abcs;
    2500   for (uInt i = 0; i < abcsd.size(); ++i) {
    2501     abcs.push_back((float)abcsd[i]);
     2519  std::vector<double> dAbcissa = getAbcissa(whichrow);
     2520  std::vector<float> abcissa;
     2521  for (uInt i = 0; i < dAbcissa.size(); ++i) {
     2522    abcissa.push_back((float)dAbcissa[i]);
    25022523  }
    25032524  std::vector<float> spec = getSpectrum(whichrow);
    25042525
    2505   fitter.setData(abcs, spec, mask);
     2526  fitter.setData(abcissa, spec, mask);
    25062527  fitter.lfit();
    25072528}
    … …  
    25662587
    25672588/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2568 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) {
     2589void 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) {
    25692590  if (outLogger || outTextFile) {
    25702591    float rms = getRms(chanMask, whichrow);
    25712592    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2593    std::vector<bool> fixed;
     2594    fixed.clear();
    25722595
    25732596    if (outLogger) {
    … …  
    25822605
    25832606/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2584 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) {
     2607void 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) {
    25852608  if (outLogger || outTextFile) {
    25862609    float rms = getRms(chanMask, whichrow);
    25872610    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2611    std::vector<bool> fixed;
     2612    fixed.clear();
    25882613
    25892614    if (outLogger) {

  • trunk/src/Scantable.h

    r2064 r2081  
    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);

  • trunk/src/ScantableWrapper.h

    r2047 r2081  
    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.