Changeset 352 for trunk/src


Ignore:
Timestamp:
02/02/05 16:31:38 (20 years ago)
Author:
vor010
Message:

SDLineFinder - further interface changes to achieve a better structure of the program

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDLineFinder.cc

    r351 r352  
    4646
    4747namespace asap {
    48 
    49 ///////////////////////////////////////////////////////////////////////////////
    50 //
    51 // IStatHolder - an abstract class to collect statistics from the running
    52 //               mean calculator, if necessary.
    53 //               We define it here, because it is used in LFRunningMean and
    54 //               SDLineFinder only
    55 //
    56 
    57 struct IStatHolder  {
    58    // This function is called for each spectral channel processed by
    59    // the running mean calculator. The order of channel numbers may be
    60    // arbitrary
    61    // ch - a number of channel, corresponding to (approximately) the centre
    62    //      of the running box
    63    // box_nchan - number of channels in the box
    64    //
    65    virtual void accumulate(int ch, Float sum, Float sum2, int box_nchan)
    66                       throw(AipsError) = 0;                   
    67 };
    68 
    69 //
    70 ///////////////////////////////////////////////////////////////////////////////
    71 
    72 ///////////////////////////////////////////////////////////////////////////////
    73 //
    74 // SignAccumulator - a simple class to deal with running mean statistics:
    75 //                   it stores the sign of the value-mean only
    76 //
    77 class SignAccumulator : public IStatHolder {
    78    Vector<Int>  sign;                // either +1, -1 or 0
    79    const Vector<Float>   &spectrum;  // a reference to the spectrum
    80                                      // to calculate the sign   
    81 public:
    82    // all channels >=nchan are ignored
    83    SignAccumulator(uInt nchan, const Vector<Float> &in_spectrum);
    84    
    85    // This function is called for each spectral channel processed by
    86    // the running mean calculator. The order of channel numbers may be
    87    // arbitrary
    88    // ch - a number of channel, corresponding to (approximately) the centre
    89    //      of the running box
    90    // box_nchan - number of channels in the box
    91    //
    92    virtual void accumulate(int ch, Float sum, Float, int box_nchan)
    93                       throw(AipsError);
    94 
    95    // access to the sign
    96    const Vector<Int>& getSigns() const throw();
    97 };
    98 
    99 //
    100 ///////////////////////////////////////////////////////////////////////////////
    10148
    10249///////////////////////////////////////////////////////////////////////////////
    … …  
    185132//                    consequtive channels. Prefix LF stands for Line Finder
    186133//
    187 class LFAboveThreshold {
     134class LFAboveThreshold : protected LFLineListOperations {
    188135   // temporary line edge channels and flag, which is True if the line
    189136   // was detected in the previous channels.
    … …  
    218165                  const casa::Vector<casa::Bool> &mask,
    219166                  const std::pair<int,int> &edge,
    220                   int max_box_nchan,
    221                   IStatHolder* statholder = NULL) throw(casa::AipsError);
     167                  int max_box_nchan) throw(casa::AipsError);
    222168
    223169protected:
    … …  
    239185
    240186} // namespace asap
    241 
    242 ///////////////////////////////////////////////////////////////////////////////
    243 //
    244 // SignAccumulator - a simple class to deal with running mean statistics:
    245 //                   it stores the sign of the value-mean only
    246 //
    247 
    248 // all channels >=nchan are ignored
    249 SignAccumulator::SignAccumulator(uInt nchan,
    250                    const Vector<Float> &in_spectrum) : sign(nchan,0),
    251                    spectrum(in_spectrum) {}
    252 
    253 
    254 // This function is called for each spectral channel processed by
    255 // the running mean calculator. The order of channel numbers may be
    256 // arbitrary
    257 // ch - a number of channel, corresponding to (approximately) the centre
    258 //      of the running box
    259 // box_nchan - number of channels in the box
    260 //
    261 void SignAccumulator::accumulate(int ch, Float sum, Float sum2, int box_nchan)
    262                    throw(AipsError)
    263 {
    264    if (ch>=sign.nelements()) return;
    265    DebugAssert(ch>=0,AipsError);
    266    DebugAssert(ch<=spectrum.nelements(), AipsError);
    267    if (box_nchan) {
    268        Float buf=spectrum[ch]-sum; ///Float(box_nchan);
    269        if (buf>0) sign[ch]=1;
    270        else if (buf<0) sign[ch]=-1;
    271        else sign[ch]=0;
    272    } else sign[ch]=0;   
    273 }
    274 
    275 // access to the sign
    276 const Vector<Int>& SignAccumulator::getSigns() const throw()
    277 {
    278    return sign;
    279 }
    280 
    281 
    282 //
    283 ///////////////////////////////////////////////////////////////////////////////
    284187
    285188///////////////////////////////////////////////////////////////////////////////
    … …  
    519422                              const casa::Vector<casa::Bool> &mask,
    520423                              const std::pair<int,int> &edge,
    521                               int max_box_nchan,
    522                               IStatHolder* statholder)
     424                              int max_box_nchan)
    523425                        throw(casa::AipsError)
    524426{
    … …  
    531433      // actual search algorithm
    532434      is_detected_before=False;
    533    
     435      Vector<Int> signs(spectrum.nelements(),0);
     436         
    534437      for (;running_box->haveMore();running_box->next()) {
    535438           const int ch=running_box->getChannel();
    536439           if (running_box->getNumberOfBoxPoints()>=minboxnchan)
    537                processChannel(mask[ch] && (fabs(spectrum[ch]-
    538                   running_box->getLinMean()) >=
     440               processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
    539441                  threshold*running_box->getLinVariance()), mask);
    540442           else processCurLine(mask); // just finish what was accumulated before
    541            if (statholder!=NULL)
    542               statholder->accumulate(running_box->getChannel(),
    543                                      running_box->getLinMean(),
    544                                      running_box->getLinVariance(),
    545                                      running_box->getNumberOfBoxPoints());
     443           const Float buf=running_box->aboveMean();
     444           if (buf>0) signs[ch]=1;
     445           else if (buf<0) signs[ch]=-1;
     446           else if (buf==0) signs[ch]=0;           
    546447      }
     448      if (lines.size())
     449          searchForWings(lines,signs,mask,edge);
    547450  }
    548451  catch (const AipsError &ae) {
    … …  
    559462///////////////////////////////////////////////////////////////////////////////
    560463//
    561 // SDLineFinder::IntersectsWith  -  An auxiliary object function to test
    562 // whether two lines have a non-void intersection
     464// LFLineListOperations::IntersectsWith  -  An auxiliary object function
     465//                to test whether two lines have a non-void intersection
    563466//
    564467
    565468
    566469// line1 - range of the first line: start channel and stop+1
    567 SDLineFinder::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
     470LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
    568471                          line1(in_line1) {}
    569472
    … …  
    572475// common channel, and false otherwise
    573476// line2 - range of the second line: start channel and stop+1
    574 bool SDLineFinder::IntersectsWith::operator()(const std::pair<int,int> &line2)
     477bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
    575478                          const throw()
    576479{
    … …  
    585488///////////////////////////////////////////////////////////////////////////////
    586489//
    587 // SDLineFinder::BuildUnion - An auxiliary object function to build a union
     490// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
    588491// of several lines to account for a possibility of merging the nearby lines
    589492//
    590493
    591494// set an initial line (can be a first line in the sequence)
    592 SDLineFinder::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
     495LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
    593496                             temp_line(line1) {}
    594497
    595498// update temp_line with a union of temp_line and new_line
    596499// provided there is no gap between the lines
    597 void SDLineFinder::BuildUnion::operator()(const std::pair<int,int> &new_line)
     500void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
    598501                                   throw()
    599502{
    … …  
    603506
    604507// return the result (temp_line)
    605 const std::pair<int,int>& SDLineFinder::BuildUnion::result() const throw()
     508const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
    606509{
    607510  return temp_line;
    … …  
    613516///////////////////////////////////////////////////////////////////////////////
    614517//
    615 // SDLineFinder::LaterThan - An auxiliary object function to test whether a
     518// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
    616519// specified line is at lower spectral channels (to preserve the order in
    617520// the line list)
    … …  
    619522
    620523// setup the line to compare with
    621 SDLineFinder::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
     524LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
    622525                         line1(in_line1) {}
    623526
    624527// return true if line2 should be placed later than line1
    625528// in the ordered list (so, it is at greater channel numbers)
    626 bool SDLineFinder::LaterThan::operator()(const std::pair<int,int> &line2)
     529bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
    627530                          const throw()
    628531{
    … …  
    742645     LFAboveThreshold lfalg(new_lines,min_nchan, threshold);
    743646
    744      SignAccumulator sacc(spectrum.nelements(),spectrum);
    745 
    746 
    747647     try {
    748          lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan,&sacc);
     648         lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
    749649     }
    750650     catch(const AipsError &ae) {
    … …  
    754654     first_pass=False;
    755655     if (!new_lines.size()) break; // nothing new
    756 
    757      searchForWings(new_lines, sacc.getSigns());
    758656
    759657     // update the list (lines) merging intervals, if necessary
    … …  
    849747}
    850748
     749//
     750///////////////////////////////////////////////////////////////////////////////
     751
     752
     753///////////////////////////////////////////////////////////////////////////////
     754//
     755// LFLineListOperations - a class incapsulating  operations with line lists
     756//                        The LF prefix stands for Line Finder
     757//
     758
    851759// concatenate two lists preserving the order. If two lines appear to
    852760// be adjacent, they are joined into the new one
    853 void SDLineFinder::addNewSearchResult(const std::list<pair<int, int> > &newlines,
     761void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
    854762                         std::list<std::pair<int, int> > &lines_list)
    855763                        throw(AipsError)
    … …  
    888796  } 
    889797  catch (const exception &ex) {
    890       throw AipsError(String("SDLineFinder::addNewSearchResult - STL error: ")+ex.what());
     798      throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
    891799  }
    892800}
    … …  
    897805// the detection threshold. If lines becomes adjacent, they are
    898806// merged together. Any masked channel stops the extension
    899 void SDLineFinder::searchForWings(std::list<std::pair<int, int> > &newlines,
    900            const casa::Vector<casa::Int> &signs) throw(casa::AipsError)
     807void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
     808           const casa::Vector<casa::Int> &signs,
     809           const casa::Vector<casa::Bool> &mask,
     810           const std::pair<int,int> &edge) throw(casa::AipsError)
    901811{
    902812  try {
    … …  
    928838  } 
    929839  catch (const exception &ex) {
    930       throw AipsError(String("SDLineFinder::extendLines - STL error: ")+ex.what());
    931   }
    932 }
     840      throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
     841  }
     842}
     843
     844//
     845///////////////////////////////////////////////////////////////////////////////

  • trunk/src/SDLineFinder.h

    r344 r352  
    5454namespace asap {
    5555
    56 // SDLineFinder  -  a class for automated spectral line search
    57 struct SDLineFinder {
    58    SDLineFinder() throw();
    59    virtual ~SDLineFinder() throw(casa::AipsError);
     56///////////////////////////////////////////////////////////////////////////////
     57//
     58// LFLineListOperations - a class incapsulating  operations with line lists
     59//                        The LF prefix stands for Line Finder
     60//
    6061
    61    // set the scan to work with (in_scan parameter), associated mask (in_mask
    62    // parameter) and the edge channel rejection (in_edge parameter)
    63    //   if in_edge has zero length, all channels chosen by mask will be used
    64    //   if in_edge has one element only, it represents the number of
    65    //      channels to drop from both sides of the spectrum
    66    //   in_edge is introduced for convinience, although all functionality
    67    //   can be achieved using a spectrum mask only   
    68    void setScan(const SDMemTableWrapper &in_scan,
    69                 const std::vector<bool> &in_mask,
    70                 const boost::python::tuple &in_edge) throw(casa::AipsError);
    71 
    72    // search for spectral lines. Number of lines found is returned
    73    int findLines() throw(casa::AipsError);
    74 
    75    // get the mask to mask out all lines that have been found (default)
    76    // if invert=true, only channels belong to lines will be unmasked
    77    // Note: all channels originally masked by the input mask (in_mask
    78    //       in setScan) or dropped out by the edge parameter (in_edge
    79    //       in setScan) are still excluded regardless on the invert option
    80    std::vector<bool> getMask(bool invert=false) const throw(casa::AipsError);
    81 
    82    // get range for all lines found. If defunits is true (default), the
    83    // same units as used in the scan will be returned (e.g. velocity
    84    // instead of channels). If defunits is false, channels will be returned
    85    std::vector<int>   getLineRanges(bool defunits=true)
    86                                 const throw(casa::AipsError);
    87 protected:
     62struct  LFLineListOperations {
    8863   // concatenate two lists preserving the order. If two lines appear to
    8964   // be adjacent or have a non-void intersection, they are joined into
    … …  
    9873   // the detection threshold. If lines becomes adjacent, they are
    9974   // merged together. Any masked channel stops the extension
    100    void searchForWings(std::list<std::pair<int, int> > &newlines,
    101                            const casa::Vector<casa::Int> &signs)
     75   static void searchForWings(std::list<std::pair<int, int> > &newlines,
     76                       const casa::Vector<casa::Int> &signs,
     77                       const casa::Vector<casa::Bool> &mask,
     78                       const std::pair<int,int> &edge)
    10279                           throw(casa::AipsError);
    103                            
     80protected:
     81           
    10482   // An auxiliary object function to test whether two lines have a non-void
    10583   // intersection
    … …  
    141119   };
    142120   
     121   
     122};
     123
     124//
     125///////////////////////////////////////////////////////////////////////////////
     126
     127///////////////////////////////////////////////////////////////////////////////
     128//
     129// SDLineFinder  -  a class for automated spectral line search
     130//
     131//
     132
     133struct SDLineFinder : protected LFLineListOperations {
     134   SDLineFinder() throw();
     135   virtual ~SDLineFinder() throw(casa::AipsError);
     136
     137   // set the scan to work with (in_scan parameter), associated mask (in_mask
     138   // parameter) and the edge channel rejection (in_edge parameter)
     139   //   if in_edge has zero length, all channels chosen by mask will be used
     140   //   if in_edge has one element only, it represents the number of
     141   //      channels to drop from both sides of the spectrum
     142   //   in_edge is introduced for convinience, although all functionality
     143   //   can be achieved using a spectrum mask only   
     144   void setScan(const SDMemTableWrapper &in_scan,
     145                const std::vector<bool> &in_mask,
     146                const boost::python::tuple &in_edge) throw(casa::AipsError);
     147
     148   // search for spectral lines. Number of lines found is returned
     149   int findLines() throw(casa::AipsError);
     150
     151   // get the mask to mask out all lines that have been found (default)
     152   // if invert=true, only channels belong to lines will be unmasked
     153   // Note: all channels originally masked by the input mask (in_mask
     154   //       in setScan) or dropped out by the edge parameter (in_edge
     155   //       in setScan) are still excluded regardless on the invert option
     156   std::vector<bool> getMask(bool invert=false) const throw(casa::AipsError);
     157
     158   // get range for all lines found. If defunits is true (default), the
     159   // same units as used in the scan will be returned (e.g. velocity
     160   // instead of channels). If defunits is false, channels will be returned
     161   std::vector<int>   getLineRanges(bool defunits=true)
     162                                const throw(casa::AipsError);
    143163private:
    144164   casa::CountedConstPtr<SDMemTable> scan; // the scan to work with
    … …  
    159179   // a buffer for the spectrum
    160180   mutable casa::Vector<casa::Float>  spectrum;
     181};
    161182
    162 };
     183//
     184///////////////////////////////////////////////////////////////////////////////
     185
    163186} // namespace asap
    164187#endif // #ifndef SDLINEFINDER_H
Note: See TracChangeset for help on using the changeset viewer.