Changeset 2852


Ignore:
Timestamp:
09/13/13 19:09:11 (11 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: Yes (CAS-5139/TRAC-290)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: A couple of private functions are added to STSideBandSep class.

Test Programs:

Put in Release Notes: Yes

Module(s): asap.sbseparator

Description: Implemented a way to transfer channel based flag of imput tables to output tables. Channel flag is accumulated by logical SUM.


Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/STSideBandSep.cpp

    r2838 r2852  
    261261  remRowIds.resize(0);
    262262  Matrix<float> specMat(nchan_, nshift_);
     263  Matrix<bool> flagMat(nchan_, nshift_);
    263264  vector<float> sigSpec(nchan_), imgSpec(nchan_);
     265  Vector<bool> flagVec(nchan_);
    264266  vector<uInt> tabIdvec;
    265267
     
    276278    // Get a set of spectra to solve
    277279    if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
    278                            specMat, tabIdvec)){
     280                           specMat, flagMat, tabIdvec)){
    279281      remRowIds.push_back(irow);
    280282#ifdef KS_DEBUG
     
    289291      // unflag the spectrum since there should be some valid data
    290292      sigTab_p->flagRow(vector<uInt>(irow), true);
     293      // need to unflag whole channels anyway
    291294      sigTab_p->flag(irow, vector<bool>(), true);
    292295    }
     296    // apply channel flag
     297    flagVec = collapseFlag(flagMat, tabIdvec, true);
     298    //boolVec = !boolVec; // flag
     299    sigTab_p->flag(irow, flagVec.tovector(), false);
    293300
    294301    // Solve image sideband
     
    299306        // unflag the spectrum since there should be some valid data
    300307        imgTab_p->flagRow(vector<uInt>(irow), true);
     308        // need to unflag whole channels anyway
    301309        imgTab_p->flag(irow, vector<bool>(), true);
    302310      }
     311      // apply channel flag
     312      flagVec = collapseFlag(flagMat, tabIdvec, false);
     313      //boolVec = !boolVec; // flag
     314      imgTab_p->flag(irow, flagVec.tovector(), true);
    303315    }
    304316  } // end of row loop
     
    594606};
    595607
     608// bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
     609//                                    const double dirX, const double dirY,
     610//                                    Matrix<float> &specMat, vector<uInt> &tabIdvec)
    596611bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
    597612                                      const double dirX, const double dirY,
    598                                       Matrix<float> &specMat, vector<uInt> &tabIdvec)
     613                                      Matrix<float> &specMat,
     614                                      Matrix<bool> &flagMat,
     615                                      vector<uInt> &tabIdvec)
    599616{
    600617  LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
     
    603620  specMat.resize(nchan_, nshift_);
    604621  Vector<float> spec;
     622  Vector<bool> boolVec;
    605623  uInt nspec = 0;
    606624  STMath stm(false); // insitu has no effect for average.
     
    675693    // Interpolate flagged channels of the spectrum.
    676694    Vector<Float> tmpSpec = avetab_p->getSpectrum(0);
     695    // Mask is true if the data is valid (!flag)
    677696    vector<bool> mask = avetab_p->getMask(0);
    678697    mathutil::doZeroOrderInterpolation(tmpSpec, mask);
    679698    spec.reference(specMat.column(nspec));
    680699    spec = tmpSpec;
     700    boolVec.reference(flagMat.column(nspec));
     701    boolVec = mask; // cast std::vector to casa::Vector
     702    boolVec = !boolVec;
    681703    tabIdvec.push_back((uInt) itab);
    682704    nspec++;
     705    //Liberate from reference
     706    spec.unique();
     707    boolVec.unique();
    683708  } // end of table loop
    684709  // Check the number of selected spectra and resize matrix.
     
    686711    //shiftSpecmat.resize(nchan_, nspec, true);
    687712    specMat.resize(nchan_, nspec, true);
     713    flagMat.resize(nchan_, nspec, true);
    688714#ifdef KS_DEBUG
    689715      cout << "Could not find corresponding rows in some tables."
     
    701727  return true;
    702728};
     729
     730
     731Vector<bool> STSideBandSep::collapseFlag(const Matrix<bool> &flagMat,
     732                                         const vector<uInt> &tabIdvec,
     733                                         const bool signal)
     734{
     735  LogIO os(LogOrigin("STSideBandSep","collapseFlag()", WHERE));
     736  if (tabIdvec.size() == 0)
     737    throw(AipsError("Internal error. Table index is not defined."));
     738  if (flagMat.ncolumn() != tabIdvec.size())
     739    throw(AipsError("Internal error. The row number of input matrix is not conformant."));
     740  if (flagMat.nrow() != nchan_)
     741    throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
     742 
     743  const size_t nspec = tabIdvec.size();
     744  vector<double> *thisShift;
     745  if (signal == otherside_) {
     746    // (solve signal && solveother = T) OR (solve image && solveother = F)
     747    thisShift = &imgShift_;
     748  } else {
     749    // (solve signal && solveother = F) OR (solve image && solveother = T)
     750    thisShift =  &sigShift_;
     751 }
     752
     753  Vector<bool> outflag(nchan_, false);
     754  double tempshift;
     755  Vector<bool> shiftvec(nchan_, false);
     756  Vector<bool> accflag(nchan_, false);
     757  uInt shiftId;
     758  for (uInt i = 0 ; i < nspec; ++i) {
     759    shiftId = tabIdvec[i];
     760    tempshift = - thisShift->at(shiftId);
     761    shiftFlag(flagMat.column(i), tempshift, shiftvec);
     762    // Now accumulate Flag
     763    for (uInt j = 0 ; j < nchan_ ; ++j)
     764      accflag[j] |= shiftvec[j];
     765  }
     766  // Shift back Flag
     767  shiftFlag(accflag, thisShift->at(0), outflag);
     768
     769  return outflag;
     770}
     771
    703772
    704773vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
     
    815884  }
    816885};
     886
     887
     888void STSideBandSep::shiftFlag(const Vector<bool> &invec,
     889                                  double shift,
     890                                  Vector<bool> &outvec)
     891{
     892  LogIO os(LogOrigin("STSideBandSep","shiftFlag()", WHERE));
     893  if (invec.size() != nchan_)
     894    throw(AipsError("Internal error. The length of input vector differs from nchan_"));
     895  if (outvec.size() != nchan_)
     896    throw(AipsError("Internal error. The length of output vector differs from nchan_"));
     897
     898#ifdef KS_DEBUG
     899  cout << "Start shifting flag for " << shift << "channels" << endl;
     900#endif
     901
     902  // shift is almost integer think it as int.
     903  // tolerance should be in 0 - 1
     904  double tolerance = 0.01;
     905  // tweak shift to be in 0 ~ nchan_-1
     906  if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
     907  if (shift < 0.) shift += nchan_;
     908  double rweight = fmod(shift, 1.);
     909  bool ruse(true), luse(true);
     910  if (rweight < 0.) rweight += 1.;
     911  if (rweight < tolerance){
     912    ruse = true;
     913    luse = false;
     914  }
     915  if (rweight > 1-tolerance){
     916    ruse = false;
     917    luse = true;
     918  }
     919  uInt lchan, rchan;
     920
     921  outvec = false;
     922  for (uInt i = 0 ; i < nchan_ ; i++) {
     923    lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
     924    if (lchan < 0.) lchan += nchan_;
     925    rchan = ( (lchan + 1) % nchan_ );
     926    outvec(lchan) |= (invec(i) && luse);
     927    outvec(rchan) |= (invec(i) && ruse);
     928  }
     929};
     930
    817931
    818932void STSideBandSep::deconvolve(Matrix<float> &specmat,
  • trunk/src/STSideBandSep.h

    r2828 r2852  
    130130                         const double refval, const double refpix,
    131131                         const double increment, const int nChan);
     132  //  bool getSpectraToSolve(const int polId, const int beamId,
     133  //                     const double dirX, const double dirY,
     134  //                     Matrix<float> &specmat, vector<uInt> &tabIdvec);
    132135  bool getSpectraToSolve(const int polId, const int beamId,
    133136                         const double dirX, const double dirY,
    134                          Matrix<float> &specmat, vector<uInt> &tabIdvec);
     137                         Matrix<float> &specMat, Matrix<bool> &flagMat,
     138                         vector<uInt> &tabIdvec);
    135139
    136   vector<float> solve(const Matrix<float> &specmat,
     140  vector<float> solve(const Matrix<float> &specMat,
    137141                      const vector<uInt> &tabIdvec,
    138142                      const bool signal = true);
    139143
     144  Vector<bool> collapseFlag(const Matrix<bool> &flagMat,
     145                            const vector<uInt> &tabIdvec,
     146                            const bool signal = true);
     147
    140148  void shiftSpectrum(const Vector<float> &invec, double shift,
    141149                     Vector<float> &outvec);
     150
     151  void shiftFlag(const Vector<bool> &invec, double shift,
     152                     Vector<bool> &outvec);
    142153
    143154  void deconvolve(Matrix<float> &specmat, const vector<double> shiftvec,
Note: See TracChangeset for help on using the changeset viewer.