Changeset 1644 for trunk/src


Ignore:
Timestamp:
10/04/09 00:53:18 (15 years ago)
Author:
Max Voronkov
Message:

line finder: added more options on how the noise is to be estimated. See doc on linefinder.set_options for more info

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/STLineFinder.cpp

    r1643 r1644  
    111111
    112112protected:
    113    // supplementary function to control running mean calculations.
    114    // It adds a specified channel to the running mean box and
     113   // supplementary function to control running mean/median calculations.
     114   // It adds a specified channel to the running box and
    115115   // removes (ch-maxboxnchan+1)'th channel from there
    116116   // Channels, for which the mask is false or index is beyond the
     
    153153                                           // last point of the detected line
    154154                                           //
     155   bool itsUseMedian;                      // true if median statistics is used
     156                                           // to determine the noise level, otherwise
     157                                           // it is the mean of the lowest 80% of deviations
     158                                           // (default)
     159   int itsNoiseSampleSize;                 // sample size used to estimate the noise statistics
     160                                           // Negative value means the whole spectrum is used (default)
    155161public:
    156162
     
    158164   LFAboveThreshold(std::list<pair<int,int> > &in_lines,
    159165                    int in_min_nchan = 3,
    160                     casa::Float in_threshold = 5) throw();
     166                    casa::Float in_threshold = 5,
     167                    bool use_median = false,
     168                    int noise_sample_size = -1) throw();
    161169   virtual ~LFAboveThreshold() throw();
    162170
     
    227235   // mean of lowest 80% of the samples
    228236   float meanLowest80Percent() const;
     237
     238   // return true if the buffer is full (i.e. statistics are representative)
     239   inline bool filledToCapacity() const { return itsBufferFull;}
    229240
    230241protected:
     
    509520}
    510521
    511 // supplementary function to control running mean calculations.
    512 // It adds a specified channel to the running mean box and
     522// supplementary function to control running mean/median calculations.
     523// It adds a specified channel to the running box and
    513524// removes (ch-max_box_nchan+1)'th channel from there
    514525// Channels, for which the mask is false or index is beyond the
     
    585596///////////////////////////////////////////////////////////////////////////////
    586597//
    587 // LFAboveThreshold - a running mean algorithm for line detection
     598// LFAboveThreshold - a running mean/median algorithm for line detection
    588599//
    589600//
     
    593604LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
    594605                                   int in_min_nchan,
    595                                    casa::Float in_threshold) throw() :
     606                                   casa::Float in_threshold,
     607                                   bool use_median,
     608                                   int noise_sample_size) throw() :
    596609             min_nchan(in_min_nchan), threshold(in_threshold),
    597              lines(in_lines), running_box(NULL) {}
     610             lines(in_lines), running_box(NULL), itsUseMedian(use_median),
     611             itsNoiseSampleSize(noise_sample_size) {}
    598612
    599613LFAboveThreshold::~LFAboveThreshold() throw()
     
    711725      // an assumption made: lines occupy a small part of the spectrum
    712726
    713       DebugAssert(edge.second-edge.first,AipsError);
    714       LFNoiseEstimator ne(edge.second-edge.first);
     727      const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) :
     728                      std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
     729      DebugAssert(noiseSampleSize,AipsError);
     730      const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
     731      LFNoiseEstimator ne(noiseSampleSize);
    715732
    716733      for (;running_box->haveMore();running_box->next()) {
    717            ne.add(running_box->getLinVariance());
    718       }
    719 
    720       const Float offline_variance = ne.meanLowest80Percent();
     734           ne.add(running_box->getLinVariance());
     735           if (ne.filledToCapacity()) {
     736               break;
     737           }
     738      }
     739
     740      Float offline_variance = -1; // just a flag that it is unset
     741           
     742      if (globalNoise) {
     743          offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
     744      }
    721745
    722746      // actual search algorithm
     
    731755                                 running_box->next()) {
    732756           const int ch=running_box->getChannel();
    733            if (running_box->getNumberOfBoxPoints()>=minboxnchan)
     757           if (!globalNoise) {
     758               // add a next point for a local noise estimate
     759               ne.add(running_box->getLinVariance());
     760           }   
     761           if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
     762               if (!globalNoise) {
     763                   offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
     764               }
     765               AlwaysAssert(offline_variance>0.,AipsError);
    734766               processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
    735767                  threshold*offline_variance), mask);
    736            else processCurLine(mask); // just finish what was accumulated before
     768           } else processCurLine(mask); // just finish what was accumulated before
    737769
    738770           signs[ch]=getAboveMeanSign();
     
    863895//              Setting a very large value doesn't usually provide
    864896//              valid detections.
    865 // in_box_size  the box size for running mean calculation. Default is
     897// in_box_size  the box size for running mean/median calculation. Default is
    866898//              1./5. of the whole spectrum size
     899// in_noise_box the box size for off-line noise estimation (if working with
     900//              local noise. Negative value means use global noise estimate
     901//              Default is -1 (i.e. estimate using the whole spectrum)
     902// in_median    true if median statistics is used as opposed to average of
     903//              the lowest 80% of deviations (default)
    867904void STLineFinder::setOptions(const casa::Float &in_threshold,
    868905                              const casa::Int &in_min_nchan,
    869906                              const casa::Int &in_avg_limit,
    870                               const casa::Float &in_box_size) throw()
     907                              const casa::Float &in_box_size,
     908                              const casa::Float &in_noise_box,
     909                              const casa::Bool &in_median) throw()
    871910{
    872911  threshold=in_threshold;
     
    874913  avg_limit=in_avg_limit;
    875914  box_size=in_box_size;
     915  itsNoiseBox = in_noise_box;
     916  itsUseMedian = in_median;
    876917}
    877918
     
    958999      throw AipsError("STLineFinder::findLines - box_size is too small");
    9591000
     1001  // number of elements in the sample for noise estimate
     1002  const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
     1003
     1004  if ((noise_box!= -1) and (noise_box<2))
     1005      throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
     1006
    9601007  spectrum.resize();
    9611008  spectrum = Vector<Float>(scan->getSpectrum(whichRow));
     
    9811028     try {
    9821029         // line find algorithm
    983          LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
     1030         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
    9841031         lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
    9851032         signs.resize(lfalg.getSigns().nelements());
  • trunk/src/STLineFinder.h

    r1353 r1644  
    150150   // in_box_size  the box size for running mean calculation. Default is
    151151   //              1./5. of the whole spectrum size
     152   // in_noise_box the box size for off-line noise estimation (if working with
     153   //              local noise. Negative value means use global noise estimate
     154   //              Default is -1 (i.e. estimate using the whole spectrum)
     155   // in_median    true if median statistics is used as opposed to average of
     156   //              the lowest 80% of deviations (default)
    152157   void setOptions(const casa::Float &in_threshold=sqrt(3.),
    153158                   const casa::Int &in_min_nchan=3,
    154159                   const casa::Int &in_avg_limit=8,
    155                    const casa::Float &in_box_size=0.2) throw();
     160                   const casa::Float &in_box_size=0.2,
     161                   const casa::Float &in_noise_box=-1.,
     162                   const casa::Bool &in_median = casa::False) throw();
    156163
    157164   // set the scan to work with (in_scan parameter)
     
    241248   // a buffer for the spectrum
    242249   mutable casa::Vector<casa::Float>  spectrum;
     250
     251   // the box size for off-line noise estimation (if working with
     252   // local noise. Negative value means use global noise estimate
     253   // Default is -1 (i.e. estimate using the whole spectrum)
     254   casa::Float itsNoiseBox;
     255
     256   // true if median statistics is used as opposed to average of
     257   // the lowest 80% of deviations (default)
     258   casa::Bool itsUseMedian;
    243259};
    244260
Note: See TracChangeset for help on using the changeset viewer.