Changeset 1644 for trunk/src

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:

• STLineFinder.cpp (modified) (13 diffs)
• STLineFinder.h (modified) (2 diffs)

trunk/src/STLineFinder.cpp

@@ -111 +111 @@
 
 protected:
-   // supplementary function to control running mean calculations.
-   // It adds a specified channel to the running mean box and
+   // supplementary function to control running mean/median calculations.
+   // It adds a specified channel to the running box and
    // removes (ch-maxboxnchan+1)'th channel from there
    // Channels, for which the mask is false or index is beyond the
@@ -153 +153 @@
                                            // last point of the detected line
                                            //
+   bool itsUseMedian;                      // true if median statistics is used
+                                           // to determine the noise level, otherwise
+                                           // it is the mean of the lowest 80% of deviations
+                                           // (default)
+   int itsNoiseSampleSize;                 // sample size used to estimate the noise statistics
+                                           // Negative value means the whole spectrum is used (default)
 public:
 
@@ -158 +164 @@
    LFAboveThreshold(std::list<pair<int,int> > &in_lines,
                     int in_min_nchan = 3,
-                    casa::Float in_threshold = 5) throw();
+                    casa::Float in_threshold = 5,
+                    bool use_median = false,
+                    int noise_sample_size = -1) throw();
    virtual ~LFAboveThreshold() throw();
 
@@ -227 +235 @@
    // mean of lowest 80% of the samples
    float meanLowest80Percent() const;
+
+   // return true if the buffer is full (i.e. statistics are representative)
+   inline bool filledToCapacity() const { return itsBufferFull;}
 
 protected:
@@ -509 +520 @@
 }
 
-// supplementary function to control running mean calculations.
-// It adds a specified channel to the running mean box and
+// supplementary function to control running mean/median calculations.
+// It adds a specified channel to the running box and
 // removes (ch-max_box_nchan+1)'th channel from there
 // Channels, for which the mask is false or index is beyond the
@@ -585 +596 @@
 ///////////////////////////////////////////////////////////////////////////////
 //
-// LFAboveThreshold - a running mean algorithm for line detection
+// LFAboveThreshold - a running mean/median algorithm for line detection
 //
 //
@@ -593 +604 @@
 LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
                                    int in_min_nchan,
-                                   casa::Float in_threshold) throw() :
+                                   casa::Float in_threshold,
+                                   bool use_median,
+                                   int noise_sample_size) throw() :
              min_nchan(in_min_nchan), threshold(in_threshold),
-             lines(in_lines), running_box(NULL) {}
+             lines(in_lines), running_box(NULL), itsUseMedian(use_median),
+             itsNoiseSampleSize(noise_sample_size) {}
 
 LFAboveThreshold::~LFAboveThreshold() throw()
@@ -711 +725 @@
       // an assumption made: lines occupy a small part of the spectrum
 
-      DebugAssert(edge.second-edge.first,AipsError);
-      LFNoiseEstimator ne(edge.second-edge.first);
+      const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) : 
+                      std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
+      DebugAssert(noiseSampleSize,AipsError);
+      const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
+      LFNoiseEstimator ne(noiseSampleSize);
 
       for (;running_box->haveMore();running_box->next()) {
-           ne.add(running_box->getLinVariance());
-      }
-
-      const Float offline_variance = ne.meanLowest80Percent();
+           ne.add(running_box->getLinVariance());
+           if (ne.filledToCapacity()) {
+               break;
+           }
+      }
+
+      Float offline_variance = -1; // just a flag that it is unset
+            
+      if (globalNoise) {
+          offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
+      }
 
       // actual search algorithm
@@ -731 +755 @@
                                  running_box->next()) {
            const int ch=running_box->getChannel();
-           if (running_box->getNumberOfBoxPoints()>=minboxnchan)
+           if (!globalNoise) {
+               // add a next point for a local noise estimate
+               ne.add(running_box->getLinVariance());
+           }    
+           if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
+               if (!globalNoise) {
+                   offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
+               }
+               AlwaysAssert(offline_variance>0.,AipsError);
                processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
                   threshold*offline_variance), mask);
-           else processCurLine(mask); // just finish what was accumulated before
+           } else processCurLine(mask); // just finish what was accumulated before
 
            signs[ch]=getAboveMeanSign();
@@ -863 +895 @@
 //              Setting a very large value doesn't usually provide
 //              valid detections.
-// in_box_size  the box size for running mean calculation. Default is
+// in_box_size  the box size for running mean/median calculation. Default is
 //              1./5. of the whole spectrum size
+// in_noise_box the box size for off-line noise estimation (if working with
+//              local noise. Negative value means use global noise estimate
+//              Default is -1 (i.e. estimate using the whole spectrum)
+// in_median    true if median statistics is used as opposed to average of
+//              the lowest 80% of deviations (default)
 void STLineFinder::setOptions(const casa::Float &in_threshold,
                               const casa::Int &in_min_nchan,
                               const casa::Int &in_avg_limit,
-                              const casa::Float &in_box_size) throw()
+                              const casa::Float &in_box_size,
+                              const casa::Float &in_noise_box,
+                              const casa::Bool &in_median) throw()
 {
   threshold=in_threshold;
@@ -874 +913 @@
   avg_limit=in_avg_limit;
   box_size=in_box_size;
+  itsNoiseBox = in_noise_box;
+  itsUseMedian = in_median;
 }
 
@@ -958 +999 @@
       throw AipsError("STLineFinder::findLines - box_size is too small");
 
+  // number of elements in the sample for noise estimate
+  const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox); 
+
+  if ((noise_box!= -1) and (noise_box<2))
+      throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
+
   spectrum.resize();
   spectrum = Vector<Float>(scan->getSpectrum(whichRow));
@@ -981 +1028 @@
      try {
          // line find algorithm
-         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
+         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
          lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
          signs.resize(lfalg.getSigns().nelements());

trunk/src/STLineFinder.h

@@ -150 +150 @@
    // in_box_size  the box size for running mean calculation. Default is
    //              1./5. of the whole spectrum size
+   // in_noise_box the box size for off-line noise estimation (if working with
+   //              local noise. Negative value means use global noise estimate
+   //              Default is -1 (i.e. estimate using the whole spectrum)
+   // in_median    true if median statistics is used as opposed to average of
+   //              the lowest 80% of deviations (default)
    void setOptions(const casa::Float &in_threshold=sqrt(3.),
                    const casa::Int &in_min_nchan=3,
                    const casa::Int &in_avg_limit=8,
-                   const casa::Float &in_box_size=0.2) throw();
+                   const casa::Float &in_box_size=0.2,
+                   const casa::Float &in_noise_box=-1.,
+                   const casa::Bool &in_median = casa::False) throw();
 
    // set the scan to work with (in_scan parameter)
@@ -241 +248 @@
    // a buffer for the spectrum
    mutable casa::Vector<casa::Float>  spectrum;
+
+   // the box size for off-line noise estimation (if working with
+   // local noise. Negative value means use global noise estimate
+   // Default is -1 (i.e. estimate using the whole spectrum)
+   casa::Float itsNoiseBox;
+
+   // true if median statistics is used as opposed to average of
+   // the lowest 80% of deviations (default)
+   casa::Bool itsUseMedian;
 };