Changeset 1644

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

Files:

• python/asaplinefind.py (modified) (3 diffs)
• src/STLineFinder.cpp (modified) (13 diffs)
• src/STLineFinder.h (modified) (2 diffs)

trunk/python/asaplinefind.py

@@ -41 +41 @@
 
     def set_options(self,threshold=1.7320508075688772,min_nchan=3,
-        avg_limit=8,box_size=0.2):
+        avg_limit=8,box_size=0.2,noise_box='all',noise_stat='mean80'):
         """
         Set the parameters of the algorithm
@@ -55 +55 @@
                           this parameter can be averaged to search for
                           broad lines. Default is 8.
-             box_size     A running mean box size specified as a fraction
+             box_size     A running mean/median box size specified as a fraction
                           of the total spectrum length. Default is 1/5
+             noise_box    Area of the spectrum used to estimate noise stats
+                          Both string values and numbers are allowed
+                          Allowed string values:
+                             'all' use all the spectrum (default)
+                             'box' noise box is the same as running mean/median
+                                   box
+                          Numeric values are defined as a fraction from the
+                          spectrum size. Values should be positive.
+                          (noise_box == box_size has the same effect as
+                           noise_box = 'box')
+             noise_stat   Statistics used to estimate noise, allowed values:
+                              'mean80' use the 80% of the lowest deviations
+                                       in the noise box (default)
+                              'median' median of deviations in the noise box
+                             
         Note:  For bad baselines threshold should be increased,
                and avg_limit decreased (or even switched off completely by
@@ -62 +77 @@
                undulations instead of real lines.
         """
-        self.finder.setoptions(threshold,min_nchan,avg_limit,box_size)
+        if noise_stat.lower() not in ["mean80",'median']:
+           raise RuntimeError, "noise_stat argument in linefinder.set_options can only be mean80 or median"
+        nStat = (noise_stat.lower() == "median")
+        nBox = -1.
+        if isinstance(noise_box,str):
+           if noise_box.lower() not in ['all','box']:
+              raise RuntimeError, "string-valued noise_box in linefinder.set_options can only be all or box"
+           if noise_box.lower() == 'box':
+              nBox = box_size
+        else:
+           nBox = float(noise_box)
+        self.finder.setoptions(threshold,min_nchan,avg_limit,box_size,nBox,nStat)
         return
 

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;
 };