Changeset 2344
- Timestamp:
- 11/08/11 17:02:32 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/Scantable.cpp
r2332 r2344 11 11 // 12 12 #include <map> 13 #include <mcheck.h> 13 14 14 15 #include <atnf/PKSIO/SrcType.h> … … 90 91 initFactories(); 91 92 setupMainTable(); 92 freqTable_ = STFrequencies(*this);93 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_ .table());93 freqTable_ = new STFrequencies(*this); 94 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_->table()); 94 95 weatherTable_ = STWeather(*this); 95 96 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table()); … … 194 195 void Scantable::copySubtables(const Scantable& other) { 195 196 Table t = table_.rwKeywordSet().asTable("FREQUENCIES"); 196 TableCopy::copyRows(t, other.freqTable_ .table());197 TableCopy::copyRows(t, other.freqTable_->table()); 197 198 t = table_.rwKeywordSet().asTable("FOCUS"); 198 199 TableCopy::copyRows(t, other.focusTable_.table()); … … 211 212 void Scantable::attachSubtables() 212 213 { 213 freqTable_ = STFrequencies(table_);214 freqTable_ = new STFrequencies(table_); 214 215 focusTable_ = STFocus(table_); 215 216 weatherTable_ = STWeather(table_); … … 222 223 Scantable::~Scantable() 223 224 { 224 //cout << "~Scantable() " << this << endl;225 delete freqTable_; 225 226 } 226 227 … … 1509 1510 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1510 1511 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1511 return freqTable_ .getSpectralCoordinate(md, mp, me, rf,1512 return freqTable_->getSpectralCoordinate(md, mp, me, rf, 1512 1513 mfreqidCol_(whichrow)); 1513 1514 } … … 1518 1519 std::vector<double> stlout; 1519 1520 int nchan = specCol_(whichrow).nelements(); 1520 String us = freqTable_ .getUnitString();1521 String us = freqTable_->getUnitString(); 1521 1522 if ( us == "" || us == "pixel" || us == "channel" ) { 1522 1523 for (int i=0; i<nchan; ++i) { … … 1585 1586 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1586 1587 SpectralCoordinate spc = 1587 freqTable_ .getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));1588 freqTable_->getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow)); 1588 1589 1589 1590 String s = "Channel"; 1590 Unit u = Unit(freqTable_ .getUnitString());1591 Unit u = Unit(freqTable_->getUnitString()); 1591 1592 if (u == Unit("km/s")) { 1592 1593 s = CoordinateUtil::axisLabel(spc, 0, True,True, True); … … 1849 1850 Double refval ; 1850 1851 Double increment ; 1851 int freqnrow = freqTable_ .table().nrow() ;1852 int freqnrow = freqTable_->table().nrow() ; 1852 1853 Vector<uInt> oldId( freqnrow ) ; 1853 1854 Vector<uInt> newId( freqnrow ) ; 1854 1855 for ( int irow = 0 ; irow < freqnrow ; irow++ ) { 1855 freqTable_ .getEntry( refpix, refval, increment, irow ) ;1856 freqTable_->getEntry( refpix, refval, increment, irow ) ; 1856 1857 /*** 1857 1858 * need to shift refpix to nmin … … 1860 1861 refval = refval - ( refpix - nmin ) * increment ; 1861 1862 refpix = 0 ; 1862 freqTable_ .setEntry( refpix, refval, increment, irow ) ;1863 freqTable_->setEntry( refpix, refval, increment, irow ) ; 1863 1864 } 1864 1865 … … 2305 2306 //fitter.setIterClipping(thresClip, nIterClip); 2306 2307 2307 int nRow = nrow();2308 std::vector<bool> chanMask;2309 2308 bool showProgress; 2310 2309 int minNRow; 2311 2310 parseProgressInfo(progressInfo, showProgress, minNRow); 2312 2311 2312 int nRow = nrow(); 2313 std::vector<bool> chanMask; 2314 2315 //-------------------------------- 2313 2316 for (int whichrow = 0; whichrow < nRow; ++whichrow) { 2314 2317 chanMask = getCompositeChanMask(whichrow, mask); 2315 2318 //fitBaseline(chanMask, whichrow, fitter); 2316 2319 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow); 2317 std::vector<int> pieceEdges ;2318 std::vector<float> params ;2320 std::vector<int> pieceEdges(nPiece+1); 2321 std::vector<float> params(nPiece*4); 2319 2322 int nClipped = 0; 2320 2323 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual); … … 2325 2328 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow); 2326 2329 } 2327 2330 //-------------------------------- 2331 2328 2332 if (outTextFile) ofs.close(); 2329 2333 … … 2394 2398 //fitBaseline(chanMask, whichrow, fitter); 2395 2399 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow); 2396 std::vector<int> pieceEdges ;2397 std::vector<float> params ;2400 std::vector<int> pieceEdges(nPiece+1); 2401 std::vector<float> params(nPiece*4); 2398 2402 int nClipped = 0; 2399 2403 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual); … … 2422 2426 2423 2427 int nChan = data.size(); 2424 std::vector<int> maskArray; 2425 std::vector<int> x; 2428 std::vector<int> maskArray(nChan); 2429 std::vector<int> x(nChan); 2430 int j = 0; 2426 2431 for (int i = 0; i < nChan; ++i) { 2427 maskArray .push_back(mask[i] ? 1 : 0);2432 maskArray[i] = mask[i] ? 1 : 0; 2428 2433 if (mask[i]) { 2429 x.push_back(i); 2430 } 2431 } 2432 2433 int initNData = x.size(); 2434 x[j] = i; 2435 j++; 2436 } 2437 } 2438 int initNData = j; 2439 2434 2440 if (initNData < nPiece) { 2435 2441 throw(AipsError("too few non-flagged channels")); … … 2437 2443 2438 2444 int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5)); 2439 std::vector<double> invEdge; 2440 idxEdge.clear(); 2441 idxEdge.push_back(x[0]); 2445 std::vector<double> invEdge(nPiece-1); 2446 idxEdge[0] = x[0]; 2442 2447 for (int i = 1; i < nPiece; ++i) { 2443 2448 int valX = x[nElement*i]; 2444 idxEdge .push_back(valX);2445 invEdge .push_back(1.0/(double)valX);2446 } 2447 idxEdge .push_back(x[x.size()-1]+1);2449 idxEdge[i] = valX; 2450 invEdge[i-1] = 1.0/(double)valX; 2451 } 2452 idxEdge[nPiece] = x[initNData-1]+1; 2448 2453 2449 2454 int nData = initNData; 2450 2455 int nDOF = nPiece + 3; //number of parameters to solve, namely, 4+(nPiece-1). 2451 2456 2452 std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual; 2457 std::vector<double> x1(nChan), x2(nChan), x3(nChan); 2458 std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan); 2459 std::vector<double> r1(nChan), residual(nChan); 2453 2460 for (int i = 0; i < nChan; ++i) { 2454 2461 double di = (double)i; 2455 2462 double dD = (double)data[i]; 2456 x1 .push_back(di);2457 x2 .push_back(di*di);2458 x3 .push_back(di*di*di);2459 z1 .push_back(dD);2460 x1z1 .push_back(dD*di);2461 x2z1 .push_back(dD*di*di);2462 x3z1 .push_back(dD*di*di*di);2463 r1 .push_back(0.0);2464 residual .push_back(0.0);2463 x1[i] = di; 2464 x2[i] = di*di; 2465 x3[i] = di*di*di; 2466 z1[i] = dD; 2467 x1z1[i] = dD*di; 2468 x2z1[i] = dD*di*di; 2469 x3z1[i] = dD*di*di*di; 2470 r1[i] = 0.0; 2471 residual[i] = 0.0; 2465 2472 } 2466 2473 … … 2540 2547 } 2541 2548 2542 std::vector<double> invDiag ;2549 std::vector<double> invDiag(nDOF); 2543 2550 for (int i = 0; i < nDOF; ++i) { 2544 invDiag .push_back(1.0/xMatrix[i][i]);2551 invDiag[i] = 1.0/xMatrix[i][i]; 2545 2552 for (int j = 0; j < nDOF; ++j) { 2546 2553 xMatrix[i][j] *= invDiag[i]; … … 2574 2581 //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of 2575 2582 //cubic terms for the other pieces (in case nPiece>1). 2576 std::vector<double> y; 2577 y.clear(); 2583 std::vector<double> y(nDOF); 2578 2584 for (int i = 0; i < nDOF; ++i) { 2579 y .push_back(0.0);2585 y[i] = 0.0; 2580 2586 for (int j = 0; j < nDOF; ++j) { 2581 2587 y[i] += xMatrix[i][nDOF+j]*zMatrix[j]; … … 2587 2593 double a2 = y[2]; 2588 2594 double a3 = y[3]; 2589 params.clear(); 2590 2595 2596 int j = 0; 2591 2597 for (int n = 0; n < nPiece; ++n) { 2592 2598 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) { 2593 2599 r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i]; 2594 residual[i] = z1[i] - r1[i]; 2595 }2596 params .push_back(a0);2597 params .push_back(a1);2598 params .push_back(a2);2599 params.push_back(a3);2600 } 2601 params[j] = a0; 2602 params[j+1] = a1; 2603 params[j+2] = a2; 2604 params[j+3] = a3; 2605 j += 4; 2600 2606 2601 2607 if (n == nPiece-1) break; … … 2607 2613 a2 += 3.0*d*iE*iE; 2608 2614 a3 -= d*iE*iE*iE; 2615 } 2616 2617 //subtract constant value for masked regions at the edge of spectrum 2618 if (idxEdge[0] > 0) { 2619 int n = idxEdge[0]; 2620 for (int i = 0; i < idxEdge[0]; ++i) { 2621 //--cubic extrapolate-- 2622 //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i]; 2623 //--linear extrapolate-- 2624 //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n]; 2625 //--constant-- 2626 r1[i] = r1[n]; 2627 } 2628 } 2629 if (idxEdge[nPiece] < nChan) { 2630 int n = idxEdge[nPiece]-1; 2631 for (int i = idxEdge[nPiece]; i < nChan; ++i) { 2632 //--cubic extrapolate-- 2633 //int m = 4*(nPiece-1); 2634 //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i]; 2635 //--linear extrapolate-- 2636 //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n]; 2637 //--constant-- 2638 r1[i] = r1[n]; 2639 } 2640 } 2641 2642 for (int i = 0; i < nChan; ++i) { 2643 residual[i] = z1[i] - r1[i]; 2609 2644 } 2610 2645 … … 2638 2673 nClipped = initNData - nData; 2639 2674 2640 std::vector<float> result ;2675 std::vector<float> result(nChan); 2641 2676 if (getResidual) { 2642 2677 for (int i = 0; i < nChan; ++i) { 2643 result .push_back((float)residual[i]);2678 result[i] = (float)residual[i]; 2644 2679 } 2645 2680 } else { 2646 2681 for (int i = 0; i < nChan; ++i) { 2647 result .push_back((float)r1[i]);2682 result[i] = (float)r1[i]; 2648 2683 } 2649 2684 } … … 2652 2687 } 2653 2688 2654 2689 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) 2655 2690 { 2656 2691 nWaves.clear();
Note:
See TracChangeset
for help on using the changeset viewer.