Changeset 881


Ignore:
Timestamp:
03/08/06 11:58:43 (18 years ago)
Author:
mar637
Message:

added linefinder

Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/SDLineFinder.cc

    r551 r881  
    11//#---------------------------------------------------------------------------
    2 //# SDLineFinder.cc: A class for automated spectral line search
     2//# STLineFinder.cc: A class for automated spectral line search
    33//#--------------------------------------------------------------------------
    44//# Copyright (C) 2004
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id:
     29//# $Id:$
    3030//#---------------------------------------------------------------------------
    3131
     
    5050///////////////////////////////////////////////////////////////////////////////
    5151//
    52 // RunningBox -    a running box calculator. This class implements 
     52// RunningBox -    a running box calculator. This class implements
    5353//                 interations over the specified spectrum and calculates
    5454//                 running box filter statistics.
     
    5757class RunningBox {
    5858   // The input data to work with. Use reference symantics to avoid
    59    // an unnecessary copying   
     59   // an unnecessary copying
    6060   const casa::Vector<casa::Float>  &spectrum; // a buffer for the spectrum
    6161   const casa::Vector<casa::Bool>   &mask; // associated mask
    6262   const std::pair<int,int>         &edge; // start and stop+1 channels
    6363                                           // to work with
    64    
     64
    6565   // statistics for running box filtering
    6666   casa::Float sumf;       // sum of fluxes
     
    6969   casa::Float sumch2;     // sum of squares of channel numbers (for linear fit)
    7070   casa::Float sumfch;     // sum of flux*(channel number) (for linear fit)
    71    
     71
    7272   int box_chan_cntr;     // actual number of channels in the box
    7373   int max_box_nchan;     // maximum allowed number of channels in the box
     
    9090                 const std::pair<int,int>         &in_edge,
    9191                 int in_max_box_nchan) throw(AipsError);
    92                  
     92
    9393   // access to the statistics
    9494   const casa::Float& getLinMean() const throw(AipsError);
     
    9696   const casa::Float aboveMean() const throw(AipsError);
    9797   int getChannel() const throw();
    98    
     98
    9999   // actual number of channels in the box (max_box_nchan, if no channels
    100100   // are masked)
     
    109109   // go to start
    110110   void rewind() throw(AipsError);
    111  
     111
    112112protected:
    113113   // supplementary function to control running mean calculations.
     
    142142                                           // the detection criterion, to be
    143143                                           // a detection
    144    casa::Float threshold;                  // detection threshold - the 
     144   casa::Float threshold;                  // detection threshold - the
    145145                                           // minimal signal to noise ratio
    146146   std::list<pair<int,int> > &lines;       // list where detections are saved
     
    157157                    casa::Float in_threshold = 5) throw();
    158158   virtual ~LFAboveThreshold() throw();
    159    
     159
    160160   // replace the detection criterion
    161161   void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
     
    198198///////////////////////////////////////////////////////////////////////////////
    199199//
    200 // RunningBox -    a running box calculator. This class implements 
     200// RunningBox -    a running box calculator. This class implements
    201201//                 interations over the specified spectrum and calculates
    202202//                 running box filter statistics.
     
    227227        ++initial_box_ch)
    228228       advanceRunningBox(initial_box_ch);
    229  
    230   if (initial_box_ch==edge.second)       
     229
     230  if (initial_box_ch==edge.second)
    231231      throw AipsError("RunningBox::rewind - too much channels are masked");
    232232
    233233  cur_channel=edge.first;
    234   start_advance=initial_box_ch-max_box_nchan/2; 
     234  start_advance=initial_box_ch-max_box_nchan/2;
    235235}
    236236
     
    291291          --box_chan_cntr;
    292292          sumf-=spectrum[ch2remove];
    293           sumf2-=square(spectrum[ch2remove]); 
     293          sumf2-=square(spectrum[ch2remove]);
    294294          sumch-=Float(ch2remove);
    295295          sumch2-=square(Float(ch2remove));
     
    319319{
    320320  AlwaysAssert(box_chan_cntr, AipsError);
    321  
     321
    322322  Float mean=sumf/Float(box_chan_cntr);
    323323
     
    385385               cur_line.second=running_box->getChannel()+1;
    386386           }
    387        } else processCurLine(mask);   
     387       } else processCurLine(mask);
    388388  }
    389389  catch (const AipsError &ae) {
    390390      throw;
    391   } 
     391  }
    392392  catch (const exception &ex) {
    393393      throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
     
    401401{
    402402  try {
    403        if (is_detected_before) {                     
     403       if (is_detected_before) {
    404404           if (cur_line.second-cur_line.first>min_nchan) {
    405405               // it was a detection. We need to change the list
    406406               Bool add_new_line=False;
    407                if (lines.size()) { 
     407               if (lines.size()) {
    408408                   for (int i=lines.back().second;i<cur_line.first;++i)
    409409                        if (mask[i]) { // one valid channel in between
     
    412412                            break;
    413413                        }
    414                } else add_new_line=True; 
    415                if (add_new_line) 
     414               } else add_new_line=True;
     415               if (add_new_line)
    416416                   lines.push_back(cur_line);
    417                else lines.back().second=cur_line.second;                   
     417               else lines.back().second=cur_line.second;
    418418           }
    419419           is_detected_before=False;
    420        }     
     420       }
    421421  }
    422422  catch (const AipsError &ae) {
    423423      throw;
    424   } 
     424  }
    425425  catch (const exception &ex) {
    426426      throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
     
    450450      running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
    451451
    452      
     452
    453453      // determine the off-line variance first
    454454      // an assumption made: lines occupy a small part of the spectrum
    455        
     455
    456456      std::vector<float> variances(edge.second-edge.first);
    457457      DebugAssert(variances.size(),AipsError);
    458      
    459       for (;running_box->haveMore();running_box->next()) 
     458
     459      for (;running_box->haveMore();running_box->next())
    460460           variances[running_box->getChannel()-edge.first]=
    461461                                running_box->getLinVariance();
    462        
     462
    463463      // in the future we probably should do a proper Chi^2 estimation
    464464      // now a simple 80% of smaller values will be used.
     
    466466      // due to a bias of the Chi^2 distribution.
    467467      stable_sort(variances.begin(),variances.end());
    468      
     468
    469469      Float offline_variance=0;
    470       uInt offline_cnt=uInt(0.8*variances.size());     
     470      uInt offline_cnt=uInt(0.8*variances.size());
    471471      if (!offline_cnt) offline_cnt=variances.size(); // no much else left,
    472472                                    // although it is very inaccurate
    473       for (uInt n=0;n<offline_cnt;++n) 
    474            offline_variance+=variances[n];         
    475       offline_variance/=Float(offline_cnt);       
    476            
     473      for (uInt n=0;n<offline_cnt;++n)
     474           offline_variance+=variances[n];
     475      offline_variance/=Float(offline_cnt);
     476
    477477      // actual search algorithm
    478478      is_detected_before=False;
     
    502502  catch (const AipsError &ae) {
    503503      throw;
    504   } 
     504  }
    505505  catch (const exception &ex) {
    506506      throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
     
    583583  if (line2.second<line1.first) return false; // line2 is at lower channels
    584584  if (line2.first>line1.second) return true; // line2 is at upper channels
    585  
     585
    586586  // line2 intersects with line1. We should have no such situation in
    587587  // practice
     
    595595///////////////////////////////////////////////////////////////////////////////
    596596//
    597 // SDLineFinder  -  a class for automated spectral line search
    598 //
    599 //
    600 
    601 SDLineFinder::SDLineFinder() throw() : edge(0,0)
     597// STLineFinder  -  a class for automated spectral line search
     598//
     599//
     600
     601STLineFinder::STLineFinder() throw() : edge(0,0)
    602602{
    603603  setOptions();
     
    614614// in_avg_limit perform the averaging of no more than in_avg_limit
    615615//              adjacent channels to search for broad lines
    616 //              Default is 8, but for a bad baseline shape this 
     616//              Default is 8, but for a bad baseline shape this
    617617//              parameter should be decreased (may be even down to a
    618618//              minimum of 1 to disable this option) to avoid
    619619//              confusing of baseline undulations with a real line.
    620 //              Setting a very large value doesn't usually provide 
    621 //              valid detections. 
     620//              Setting a very large value doesn't usually provide
     621//              valid detections.
    622622// in_box_size  the box size for running mean calculation. Default is
    623623//              1./5. of the whole spectrum size
    624 void SDLineFinder::setOptions(const casa::Float &in_threshold,
     624void STLineFinder::setOptions(const casa::Float &in_threshold,
    625625                              const casa::Int &in_min_nchan,
    626626                              const casa::Int &in_avg_limit,
     
    633633}
    634634
    635 SDLineFinder::~SDLineFinder() throw(AipsError) {}
     635STLineFinder::~STLineFinder() throw(AipsError) {}
    636636
    637637// set scan to work with (in_scan parameter), associated mask (in_mask
     
    641641//      channels to drop from both sides of the spectrum
    642642//   in_edge is introduced for convinience, although all functionality
    643 //   can be achieved using a spectrum mask only   
    644 void SDLineFinder::setScan(const SDMemTableWrapper &in_scan,
     643//   can be achieved using a spectrum mask only
     644void STLineFinder::setScan(const ScantableWrapper &in_scan,
    645645               const std::vector<bool> &in_mask,
    646646               const boost::python::tuple &in_edge) throw(AipsError)
     
    651651
    652652       mask=in_mask;
    653        if (mask.nelements()!=scan->nChan())
    654            throw AipsError("SDLineFinder::setScan - in_scan and in_mask have different"
     653       if (mask.nelements()!=scan->nchan())
     654           throw AipsError("STLineFinder::setScan - in_scan and in_mask have different"
    655655                           "number of spectral channels.");
    656656
     
    658658       int n=extract<int>(in_edge.attr("__len__")());
    659659       if (n>2 || n<0)
    660            throw AipsError("SDLineFinder::setScan - the length of the in_edge parameter"
     660           throw AipsError("STLineFinder::setScan - the length of the in_edge parameter"
    661661                           "should not exceed 2");
    662662       if (!n) {
    663            // all spectrum, no rejection
     663           // all spectra, no rejection
    664664           edge.first=0;
    665            edge.second=scan->nChan();
     665           edge.second=scan->nchan();
    666666       } else {
    667667           edge.first=extract<int>(in_edge.attr("__getitem__")(0));
    668668           if (edge.first<0)
    669                throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
     669               throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
    670670                                "number of channels to drop");
    671            if (edge.first>=scan->nChan())
    672                throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
     671           if (edge.first>=scan->nchan())
     672               throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
    673673           if (n==2) {
    674674               edge.second=extract<int>(in_edge.attr("__getitem__")(1));
    675675               if (edge.second<0)
    676                    throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
     676                   throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
    677677                                   "number of channels to drop");
    678                edge.second=scan->nChan()-edge.second;
    679            } else edge.second=scan->nChan()-edge.first;
     678               edge.second=scan->nchan()-edge.second;
     679           } else edge.second=scan->nchan()-edge.first;
    680680           if (edge.second<0 || (edge.first>=edge.second))
    681                throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
    682        }       
     681               throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
     682       }
    683683  }
    684684  catch(const AipsError &ae) {
    685685       // setScan is unsuccessfull, reset scan/mask/edge
    686        scan=CountedConstPtr<SDMemTable>(); // null pointer
     686       scan=CountedConstPtr<Scantable>(); // null pointer
    687687       mask.resize(0);
    688688       edge=pair<int,int>(0,0);
     
    692692
    693693// search for spectral lines. Number of lines found is returned
    694 int SDLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
    695 { 
     694int STLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
     695{
    696696  const int minboxnchan=4;
    697697  if (scan.null())
    698       throw AipsError("SDLineFinder::findLines - a scan should be set first,"
     698      throw AipsError("STLineFinder::findLines - a scan should be set first,"
    699699                      " use set_scan");
    700   DebugAssert(mask.nelements()==scan->nChan(), AipsError);
    701   int max_box_nchan=int(scan->nChan()*box_size); // number of channels in running
     700  DebugAssert(mask.nelements()==scan->nchan(), AipsError);
     701  int max_box_nchan=int(scan->nchan()*box_size); // number of channels in running
    702702                                                 // box
    703703  if (max_box_nchan<2)
    704       throw AipsError("SDLineFinder::findLines - box_size is too small");
    705 
    706   scan->getSpectrum(spectrum, whichRow);
     704      throw AipsError("STLineFinder::findLines - box_size is too small");
     705
     706  spectrum.resize();
     707  spectrum = Vector<Float>(scan->getSpectrum(whichRow));
    707708
    708709  lines.resize(0); // search from the scratch
     
    722723  while (true) {
    723724     // a buffer for new lines found at this iteration
    724      std::list<pair<int,int> > new_lines;     
     725     std::list<pair<int,int> > new_lines;
    725726
    726727     try {
     
    745746         subtractBaseline(temp_mask,9);
    746747         averageAdjacentChannels(temp_mask,avg_factor);
    747          continue; 
     748         continue;
    748749     }
    749750     keepStrongestOnly(temp_mask,new_lines,max_box_nchan);
     
    751752     addNewSearchResult(new_lines,lines);
    752753     // get a new mask
    753      temp_mask=getMask();     
    754   }
    755  
     754     temp_mask=getMask();
     755  }
     756
    756757  // an additional search for wings because in the presence of very strong
    757758  // lines temporary mean used at each iteration will be higher than
    758759  // the true mean
    759  
     760
    760761  if (lines.size())
    761762      LFLineListOperations::searchForWings(lines,signs,mask,edge);
    762      
     763
    763764  return int(lines.size());
    764765}
     
    767768// spectrum. It uses the SDFitter class. This action is required before
    768769// reducing the spectral resolution if the baseline shape is bad
    769 void SDLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
     770void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
    770771                      const casa::Int &order) throw(casa::AipsError)
    771772{
     
    783784  sdf.setExpression("poly",order);
    784785  if (!sdf.fit()) return; // fit failed, use old spectrum
    785   spectrum=casa::Vector<casa::Float>(sdf.getResidual());   
     786  spectrum=casa::Vector<casa::Float>(sdf.getResidual());
    786787}
    787788
     
    792793// array specified, rather than the field of this class
    793794// boxsize - a number of adjacent channels to average
    794 void SDLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
     795void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
    795796                                   const casa::Int &boxsize)
    796797                            throw(casa::AipsError)
     
    798799  DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError);
    799800  DebugAssert(boxsize!=0,AipsError);
    800  
     801
    801802  for (int n=edge.first;n<edge.second;n+=boxsize) {
    802803       DebugAssert(n<spectrum.nelements(),AipsError);
     
    807808                mean+=spectrum[k];
    808809                ++nboxch;
    809             }     
     810            }
    810811       if (nboxch<boxsize) // mask these channels
    811812           for (int k=n;k<n+boxsize && k<edge.second;++k)
     
    825826//       in setScan) or dropped out by the edge parameter (in_edge
    826827//       in setScan) are still excluded regardless on the invert option
    827 std::vector<bool> SDLineFinder::getMask(bool invert)
     828std::vector<bool> STLineFinder::getMask(bool invert)
    828829                                        const throw(casa::AipsError)
    829830{
    830831  try {
    831832       if (scan.null())
    832            throw AipsError("SDLineFinder::getMask - a scan should be set first,"
     833           throw AipsError("STLineFinder::getMask - a scan should be set first,"
    833834                      " use set_scan followed by find_lines");
    834        DebugAssert(mask.nelements()==scan->nChan(), AipsError);
     835       DebugAssert(mask.nelements()==scan->nchan(), AipsError);
    835836       /*
    836837       if (!lines.size())
    837            throw AipsError("SDLineFinder::getMask - one have to search for "
     838           throw AipsError("STLineFinder::getMask - one have to search for "
    838839                           "lines first, use find_lines");
    839        */                         
     840       */
    840841       std::vector<bool> res_mask(mask.nelements());
    841842       // iterator through lines
    842843       std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
    843        for (int ch=0;ch<res_mask.size();++ch) 
     844       for (int ch=0;ch<res_mask.size();++ch)
    844845            if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
    845846            else if (!mask[ch]) res_mask[ch]=false;
    846             else {           
     847            else {
    847848                    res_mask[ch]=!invert; // no line by default
    848849                    if (cli==lines.end()) continue;
     
    852853                        ++cli; // next line in the list
    853854                 }
    854        
     855
    855856       return res_mask;
    856857  }
    857858  catch (const AipsError &ae) {
    858859      throw;
    859   } 
     860  }
    860861  catch (const exception &ex) {
    861       throw AipsError(String("SDLineFinder::getMask - STL error: ")+ex.what());
     862      throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
    862863  }
    863864}
     
    865866// get range for all lines found. The same units as used in the scan
    866867// will be returned (e.g. velocity instead of channels).
    867 std::vector<double> SDLineFinder::getLineRanges()
     868std::vector<double> STLineFinder::getLineRanges()
    868869                             const throw(casa::AipsError)
    869870{
     
    877878  for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
    878879       if (uInt(*cri)>=vel.size())
    879           throw AipsError("SDLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
     880          throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
    880881       else *outi=vel[*cri];
    881882  return res;
     
    884885// The same as getLineRanges, but channels are always used to specify
    885886// the range
    886 std::vector<int> SDLineFinder::getLineRangesInChannels()
     887std::vector<int> STLineFinder::getLineRangesInChannels()
    887888                                   const throw(casa::AipsError)
    888889{
    889890  try {
    890891       if (scan.null())
    891            throw AipsError("SDLineFinder::getLineRangesInChannels - a scan should be set first,"
     892           throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
    892893                      " use set_scan followed by find_lines");
    893        DebugAssert(mask.nelements()==scan->nChan(), AipsError);
    894        
     894       DebugAssert(mask.nelements()==scan->nchan(), AipsError);
     895
    895896       if (!lines.size())
    896            throw AipsError("SDLineFinder::getLineRangesInChannels - one have to search for "
     897           throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
    897898                           "lines first, use find_lines");
    898                            
     899
    899900       std::vector<int> res(2*lines.size());
    900901       // iterator through lines & result
    901902       std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
    902903       std::vector<int>::iterator ri=res.begin();
    903        for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {           
     904       for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
    904905            *ri=cli->first;
    905             if (++ri!=res.end()) 
    906                 *ri=cli->second-1;         
    907        }       
     906            if (++ri!=res.end())
     907                *ri=cli->second-1;
     908       }
    908909       return res;
    909910  }
    910911  catch (const AipsError &ae) {
    911912      throw;
    912   } 
     913  }
    913914  catch (const exception &ex) {
    914       throw AipsError(String("SDLineFinder::getLineRanges - STL error: ")+ex.what());
     915      throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what());
    915916  }
    916917}
     
    920921// an auxiliary function to remove all lines from the list, except the
    921922// strongest one (by absolute value). If the lines removed are real,
    922 // they will be find again at the next iteration. This approach 
    923 // increases the number of iterations required, but is able to remove 
     923// they will be find again at the next iteration. This approach
     924// increases the number of iterations required, but is able to remove
    924925// the sidelobes likely to occur near strong lines.
    925926// Later a better criterion may be implemented, e.g.
    926927// taking into consideration the brightness of different lines. Now
    927 // use the simplest solution     
     928// use the simplest solution
    928929// temp_mask - mask to work with (may be different from original mask as
    929930// the lines previously found may be masked)
     
    931932//                 nothing will be done if it is empty
    932933// max_box_nchan - channels in the running box for baseline filtering
    933 void SDLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
     934void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
    934935                  std::list<std::pair<int, int> > &lines2update,
    935936                  int max_box_nchan)
     
    952953
    953954           if (li==lines2update.end()) break; // no more lines
    954            const int ch=running_box.getChannel();         
     955           const int ch=running_box.getChannel();
    955956           if (ch>=li->first && ch<li->second)
    956957               if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean()))
     
    964965               tmp_flux=-1;
    965966           }
    966       }     
     967      }
    967968      std::list<std::pair<int,int> > res;
    968969      res.splice(res.end(),lines2update,strongli);
     
    972973  catch (const AipsError &ae) {
    973974      throw;
    974   } 
     975  }
    975976  catch (const exception &ex) {
    976       throw AipsError(String("SDLineFinder::keepStrongestOnly - STL error: ")+ex.what());
     977      throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what());
    977978  }
    978979
     
    992993// be adjacent, they are joined into the new one
    993994void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
    994                          std::list<std::pair<int, int> > &lines_list) 
     995                         std::list<std::pair<int, int> > &lines_list)
    995996                        throw(AipsError)
    996997{
     
    998999      for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
    9991000           cli!=newlines.end();++cli) {
    1000            
     1001
    10011002           // the first item, which has a non-void intersection or touches
    10021003           // the new line
    10031004           std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
    1004                           lines_list.end(), IntersectsWith(*cli));           
    1005            // the last such item         
     1005                          lines_list.end(), IntersectsWith(*cli));
     1006           // the last such item
    10061007           std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
    10071008                          lines_list.end(), not1(IntersectsWith(*cli)));
     
    10141015           lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
    10151016
    1016            // build a union of all intersecting lines 
     1017           // build a union of all intersecting lines
    10171018           pair<int,int> union_line=for_each(lines_buffer.begin(),
    10181019                   lines_buffer.end(),BuildUnion(*cli)).result();
    1019            
     1020
    10201021           // search for a right place for the new line (union_line) and add
    10211022           std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
     
    10261027  catch (const AipsError &ae) {
    10271028      throw;
    1028   } 
     1029  }
    10291030  catch (const exception &ex) {
    10301031      throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
     
    10501051                if (signs[n]==signs[li->first] && signs[li->first])
    10511052                    li->first=n;
    1052                 else break;   
     1053                else break;
    10531054           }
    10541055           // update the right hand side
     
    10571058                if (signs[n]==signs[li->second-1] && signs[li->second-1])
    10581059                    li->second=n;
    1059                 else break;   
     1060                else break;
    10601061           }
    10611062      }
     
    10681069  catch (const AipsError &ae) {
    10691070      throw;
    1070   } 
     1071  }
    10711072  catch (const exception &ex) {
    10721073      throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
  • trunk/src/SDLineFinder.h

    r370 r881  
    11//#---------------------------------------------------------------------------
    2 //# SDLineFinder.h: A class for automated spectral line search
     2//# STLineFinder.h: A class for automated spectral line search
    33//#---------------------------------------------------------------------------
    44//# Copyright (C) 2004
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id:
     29//# $Id:$
    3030//#---------------------------------------------------------------------------
    31 #ifndef SDLINEFINDER_H
    32 #define SDLINEFINDER_H
     31#ifndef STLINEFINDER_H
     32#define STLINEFINDER_H
    3333
    3434// STL
     
    4949
    5050// ASAP
    51 #include "SDMemTableWrapper.h"
    52 #include "SDMemTable.h"
     51#include "ScantableWrapper.h"
     52#include "Scantable.h"
    5353
    5454namespace asap {
     
    6060//
    6161
    62 struct  LFLineListOperations {
     62struct LFLineListOperations {
    6363   // concatenate two lists preserving the order. If two lines appear to
    64    // be adjacent or have a non-void intersection, they are joined into 
     64   // be adjacent or have a non-void intersection, they are joined into
    6565   // the new line
    6666   static void addNewSearchResult(const std::list<std::pair<int, int> >
     
    7979                           throw(casa::AipsError);
    8080protected:
    81            
     81
    8282   // An auxiliary object function to test whether two lines have a non-void
    8383   // intersection
     
    105105        const std::pair<int,int>& result() const throw();
    106106   };
    107    
     107
    108108   // An auxiliary object function to test whether a specified line
    109109   // is at lower spectral channels (to preserve the order in the line list)
     
    117117        // in the ordered list (so, it is at greater channel numbers)
    118118        bool operator()(const std::pair<int,int> &line2) const throw();
    119    }; 
    120    
    121    
     119   };
     120
     121
    122122};
    123123
     
    127127///////////////////////////////////////////////////////////////////////////////
    128128//
    129 // SDLineFinder  -  a class for automated spectral line search
    130 //
    131 //
    132 
    133 struct SDLineFinder : protected LFLineListOperations {
    134    SDLineFinder() throw();
    135    virtual ~SDLineFinder() throw(casa::AipsError);
     129// STLineFinder  -  a class for automated spectral line search
     130//
     131//
     132
     133struct STLineFinder : protected LFLineListOperations {
     134   STLineFinder() throw();
     135   virtual ~STLineFinder() throw(casa::AipsError);
    136136
    137137   // set the parameters controlling algorithm
     
    145145   // in_avg_limit perform the averaging of no more than in_avg_limit
    146146   //              adjacent channels to search for broad lines
    147    //              Default is 8, but for a bad baseline shape this 
     147   //              Default is 8, but for a bad baseline shape this
    148148   //              parameter should be decreased (may be even down to a
    149149   //              minimum of 1 to disable this option) to avoid
    150150   //              confusing of baseline undulations with a real line.
    151    //              Setting a very large value doesn't usually provide 
    152    //              valid detections. 
     151   //              Setting a very large value doesn't usually provide
     152   //              valid detections.
    153153   // in_box_size  the box size for running mean calculation. Default is
    154154   //              1./5. of the whole spectrum size
     
    164164   //      channels to drop from both sides of the spectrum
    165165   //   in_edge is introduced for convinience, although all functionality
    166    //   can be achieved using a spectrum mask only   
    167    void setScan(const SDMemTableWrapper &in_scan,
     166   //   can be achieved using a spectrum mask only
     167   void setScan(const ScantableWrapper &in_scan,
    168168                const std::vector<bool> &in_mask,
    169169                const boost::python::tuple &in_edge) throw(casa::AipsError);
     
    171171   // search for spectral lines for a row specified by whichRow and
    172172   // Beam/IF/Pol specified by current cursor set for the scantable
    173    // Number of lines found is returned   
     173   // Number of lines found is returned
    174174   int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError);
    175175
     
    182182
    183183   // get range for all lines found. The same units as used in the scan
    184    // will be returned (e.g. velocity instead of channels).   
     184   // will be returned (e.g. velocity instead of channels).
    185185   std::vector<double>   getLineRanges() const throw(casa::AipsError);
    186186   // The same as getLineRanges, but channels are always used to specify
     
    203203   void subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
    204204                         const casa::Int &order) throw(casa::AipsError);
    205    
     205
    206206   // an auxiliary function to remove all lines from the list, except the
    207207   // strongest one (by absolute value). If the lines removed are real,
    208    // they will be find again at the next iteration. This approach 
    209    // increases the number of iterations required, but is able to remove 
     208   // they will be find again at the next iteration. This approach
     209   // increases the number of iterations required, but is able to remove
    210210   // the sidelobes likely to occur near strong lines.
    211211   // Later a better criterion may be implemented, e.g.
    212212   // taking into consideration the brightness of different lines. Now
    213    // use the simplest solution     
     213   // use the simplest solution
    214214   // temp_mask - mask to work with (may be different from original mask as
    215215   // the lines previously found may be masked)
     
    222222                                      throw (casa::AipsError);
    223223private:
    224    casa::CountedConstPtr<SDMemTable> scan; // the scan to work with
     224   casa::CountedConstPtr<Scantable> scan; // the scan to work with
    225225   casa::Vector<casa::Bool> mask;          // associated mask
    226226   std::pair<int,int> edge;                // start and stop+1 channels
    227227                                           // to work with
    228    casa::Float threshold;                  // detection threshold - the 
     228   casa::Float threshold;                  // detection threshold - the
    229229                                           // minimal signal to noise ratio
    230230   casa::Double box_size;                  // size of the box for running
     
    251251
    252252} // namespace asap
    253 #endif // #ifndef SDLINEFINDER_H
     253#endif // #ifndef STLINEFINDER_H
  • trunk/src/python_SDLineFinder.cc

    r370 r881  
    11//#---------------------------------------------------------------------------
    2 //# python_SDLineFinder.cc: python exposure of C++ SDLineFinder class
     2//# python_STLineFinder.cc: python exposure of C++ STLineFinder class
    33//#---------------------------------------------------------------------------
    44//# Copyright (C) 2004
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id:
     29//# $Id:$
    3030//#---------------------------------------------------------------------------
    3131#include <boost/python.hpp>
     
    3737namespace asap {
    3838  namespace python {
    39      void python_SDLineFinder() {
    40        class_<SDLineFinder>("linefinder")
     39     void python_STLineFinder() {
     40       class_<STLineFinder>("linefinder")
    4141         .def( init <> () )
    42          .def("setoptions",&SDLineFinder::setOptions)
    43          .def("setscan",&SDLineFinder::setScan)
    44          .def("findlines",&SDLineFinder::findLines)
    45          .def("getmask",&SDLineFinder::getMask)
    46          .def("getlineranges",&SDLineFinder::getLineRanges)
    47          .def("getlinerangesinchannels",&SDLineFinder::getLineRangesInChannels)
     42         .def("setoptions",&STLineFinder::setOptions)
     43         .def("setscan",&STLineFinder::setScan)
     44         .def("findlines",&STLineFinder::findLines)
     45         .def("getmask",&STLineFinder::getMask)
     46         .def("getlineranges",&STLineFinder::getLineRanges)
     47         .def("getlinerangesinchannels",&STLineFinder::getLineRangesInChannels)
    4848       ;
    4949     };
  • trunk/src/python_asap.cpp

    r875 r881  
    6161  asap::python::python_STMath();
    6262  asap::python::python_SDFitter();
     63  asap::python::python_STLineFinder();
    6364  /*
    6465  asap::python::python_SDWriter();
    6566  asap::python::python_SDFitTable();
    66   asap::python::python_SDLineFinder();
    6767  */
    6868  asap::python::python_SDLog();
  • trunk/src/python_asap.h

    r875 r881  
    4242    void python_STMath();
    4343    void python_SDFitter();
     44    void python_STLineFinder();
    4445    /*
    4546    void python_SDWriter();
    46     void python_SDMath();
    4747    void python_SDFitTable();
    48     void python_SDLineFinder();
    4948    */
    5049    void python_SDLog();
Note: See TracChangeset for help on using the changeset viewer.