Changeset 344 for trunk/src

Timestamp:

02/01/05 16:17:25 (20 years ago)

Author:

vor010

Message:

An algorithm to detect line wings, which are
below detection threshold, has been added

Location:

trunk/src

Files:

• SDLineFinder.cc (modified) (13 diffs)
• SDLineFinder.h (modified) (1 diff)

trunk/src/SDLineFinder.cc

@@ -44 +44 @@
 using namespace boost::python;
 
+namespace asap {
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// IStatHolder - an abstract class to collect statistics from the running
+//               mean calculator, if necessary.
+//               We define it here, because it is used in LFRunningMean and
+//               SDLineFinder only
+//
+
+struct IStatHolder  {
+   // This function is called for each spectral channel processed by
+   // the running mean calculator. The order of channel numbers may be
+   // arbitrary
+   // ch - a number of channel, corresponding to (approximately) the centre
+   //      of the running box
+   // box_nchan - number of channels in the box
+   //
+   virtual void accumulate(int ch, Float sum, Float sum2, int box_nchan)
+                      throw(AipsError) = 0;                   
+};
+
+//
+///////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// SignAccumulator - a simple class to deal with running mean statistics:
+//                   it stores the sign of the value-mean only
+//
+class SignAccumulator : public IStatHolder {
+   Vector<Int>  sign;                // either +1, -1 or 0
+   const Vector<Float>   &spectrum;  // a reference to the spectrum
+                                     // to calculate the sign   
+public:
+   // all channels >=nchan are ignored
+   SignAccumulator(uInt nchan, const Vector<Float> &in_spectrum);
+   
+   // This function is called for each spectral channel processed by
+   // the running mean calculator. The order of channel numbers may be
+   // arbitrary
+   // ch - a number of channel, corresponding to (approximately) the centre
+   //      of the running box
+   // box_nchan - number of channels in the box
+   //
+   virtual void accumulate(int ch, Float sum, Float, int box_nchan)
+                      throw(AipsError);
+
+   // access to the sign
+   const Vector<Int>& getSigns() const throw();
+};
+
+//
+///////////////////////////////////////////////////////////////////////////////
+
 ///////////////////////////////////////////////////////////////////////////////
 //
@@ -49 +104 @@
 //
 //
-
-namespace asap {
 
 // An auxiliary class implementing one pass of the line search algorithm,
@@ -96 +149 @@
 
    // find spectral lines and add them into list
-   void findLines(std::list<pair<int,int> > &lines) throw(casa::AipsError);
+   // if statholder is not NULL, the accumulate function of it will be
+   // called for each channel to save statistics
+   void findLines(std::list<pair<int,int> > &lines,
+                  IStatHolder* statholder = NULL) throw(casa::AipsError);
    
 protected:
@@ -123 +179 @@
 
 };
+
+//
+///////////////////////////////////////////////////////////////////////////////
 } // namespace asap
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// SignAccumulator - a simple class to deal with running mean statistics:
+//                   it stores the sign of the value-mean only
+//
+
+// all channels >=nchan are ignored
+SignAccumulator::SignAccumulator(uInt nchan,
+                   const Vector<Float> &in_spectrum) : sign(nchan,0),
+                   spectrum(in_spectrum) {}
+
+
+// This function is called for each spectral channel processed by
+// the running mean calculator. The order of channel numbers may be
+// arbitrary
+// ch - a number of channel, corresponding to (approximately) the centre
+//      of the running box
+// box_nchan - number of channels in the box
+//
+void SignAccumulator::accumulate(int ch, Float sum, Float, int box_nchan)
+                   throw(AipsError)
+{
+   if (ch>=sign.nelements()) return;
+   DebugAssert(ch>=0,AipsError);
+   DebugAssert(ch<=spectrum.nelements(), AipsError);
+   if (box_nchan) {
+       Float buf=spectrum[ch]-sum/Float(box_nchan);
+       if (buf>0) sign[ch]=1;
+       else if (buf<0) sign[ch]=-1;
+       else sign[ch]=0;
+   } else sign[ch]=0;    
+}
+
+// access to the sign
+const Vector<Int>& SignAccumulator::getSigns() const throw()
+{
+   return sign;
+}
+
 
 //
@@ -253 +352 @@
 
 // find spectral lines and add them into list
-void LFRunningMean::findLines(std::list<pair<int,int> > &lines)
+void LFRunningMean::findLines(std::list<pair<int,int> > &lines,
+                              IStatHolder* statholder)
                         throw(casa::AipsError)
 {
@@ -273 +373 @@
   is_detected_before=False;
 
+  if (statholder!=NULL)
+      for (int n=0;n<initial_box_ch-max_box_nchan/2;++n)
+           statholder->accumulate(n,sum,sumsq,box_chan_cntr);
+
   if (box_chan_cntr>=minboxnchan) 
       // there is a minimum amount of data. We can search in the
@@ -283 +387 @@
   for (int n=initial_box_ch-max_box_nchan/2;n<edge.second;++n) {
       advanceRunningBox(n+max_box_nchan/2); // update running mean & variance
+      if (statholder!=NULL)
+          statholder->accumulate(n,sum,sumsq,box_chan_cntr);
       if (box_chan_cntr>=minboxnchan) // have enough data to process
           processChannel(lines,n);
       else processCurLine(lines); // just finish what was accumulated before
-  }  
+  }
+  if (statholder!=NULL)
+      for (int n=edge.second;n<spectrum.nelements();++n)
+           statholder->accumulate(n,sum,sumsq,box_chan_cntr);
 }
 
@@ -468 +577 @@
   lines.resize(0); // search from the scratch
   Vector<Bool> temp_mask(mask);
-  Bool need2iterate=True;
-  while (need2iterate) {
+  
+  while (true) {
      // line find algorithm
      LFRunningMean lfalg(spectrum,temp_mask,edge,max_box_nchan,min_nchan,
                          threshold);
+     SignAccumulator sacc(spectrum.nelements(),spectrum);
+
+     // a buffer for new lines found at this iteration
      std::list<pair<int,int> > new_lines;
-     lfalg.findLines(new_lines);
-     if (!new_lines.size()) need2iterate=False;
-     temp_mask=getMask();
+     
+     lfalg.findLines(new_lines,&sacc);
+     if (!new_lines.size()) break; // nothing new
+
+     searchForWings(new_lines, sacc.getSigns());
+
      // update the list (lines) merging intervals, if necessary
-     addNewSearchResult(new_lines);
+     addNewSearchResult(new_lines,lines);
+     // get a new mask
+     temp_mask=getMask();     
   }
   return int(lines.size());
@@ -569 +686 @@
 // concatenate two lists preserving the order. If two lines appear to
 // be adjacent, they are joined into the new one
-void SDLineFinder::addNewSearchResult(const std::list<pair<int, int> > &newlines)
+void SDLineFinder::addNewSearchResult(const std::list<pair<int, int> > &newlines,
+                         std::list<std::pair<int, int> > &lines_list) 
                         throw(AipsError)
 {
@@ -578 +696 @@
            // the first item, which has a non-void intersection or touches
            // the new line
-           std::list<pair<int,int> >::iterator pos_beg=find_if(lines.begin(),
-                          lines.end(), IntersectsWith(*cli));           
+           std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
+                          lines_list.end(), IntersectsWith(*cli));           
            // the last such item          
            std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
-                          lines.end(), not1(IntersectsWith(*cli)));
+                          lines_list.end(), not1(IntersectsWith(*cli)));
 
            // extract all lines which intersect or touch a new one into
@@ -589 +707 @@
            // one.
            std::list<pair<int,int> > lines_buffer;
-           lines_buffer.splice(lines_buffer.end(),lines, pos_beg, pos_end);
+           lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
 
            // build a union of all intersecting lines 
@@ -596 +714 @@
            
            // search for a right place for the new line (union_line) and add
-           std::list<pair<int,int> >::iterator pos2insert=find_if(lines.begin(),
-                          lines.end(), LaterThan(union_line));
-           lines.insert(pos2insert,union_line);
+           std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
+                          lines_list.end(), LaterThan(union_line));
+           lines_list.insert(pos2insert,union_line);
       }
   }
@@ -608 +726 @@
   }
 }
+
+// extend all line ranges to the point where a value stored in the
+// specified vector changes (e.g. value-mean change its sign)
+// This operation is necessary to include line wings, which are below
+// the detection threshold. If lines becomes adjacent, they are
+// merged together. Any masked channel stops the extension
+void SDLineFinder::searchForWings(std::list<std::pair<int, int> > &newlines,
+           const casa::Vector<casa::Int> &signs) throw(casa::AipsError)
+{
+  try {
+      for (std::list<pair<int,int> >::iterator li=newlines.begin();
+           li!=newlines.end();++li) {
+           // update the left hand side
+           for (int n=li->first-1;n>=edge.first;--n) {
+                if (!mask[n]) break;
+                if (signs[n]==signs[li->first] && signs[li->first])
+                    li->first=n;
+                else break;    
+           }
+           // update the right hand side
+           for (int n=li->second;n<edge.second;++n) {
+                if (!mask[n]) break;
+                if (signs[n]==signs[li->second-1] && signs[li->second-1])
+                    li->second=n;
+                else break;    
+           }
+      }
+      // need to search for possible mergers.
+      std::list<std::pair<int, int> >  result_buffer;
+      addNewSearchResult(newlines,result_buffer);
+      newlines.clear();
+      newlines.splice(newlines.end(),result_buffer);
+  }
+  catch (const AipsError &ae) {
+      throw;
+  }  
+  catch (const exception &ex) {
+      throw AipsError(String("SDLineFinder::extendLines - STL error: ")+ex.what());
+  }
+}

trunk/src/SDLineFinder.h

@@ -89 +89 @@
    // be adjacent or have a non-void intersection, they are joined into 
    // the new line
-   void addNewSearchResult(const std::list<std::pair<int, int> > &newlines)
+   static void addNewSearchResult(const std::list<std::pair<int, int> >
+                  &newlines, std::list<std::pair<int, int> > &lines_list)
                            throw(casa::AipsError);
+
+   // extend all line ranges to the point where a value stored in the
+   // specified vector changes (e.g. value-mean change its sign)
+   // This operation is necessary to include line wings, which are below
+   // the detection threshold. If lines becomes adjacent, they are
+   // merged together. Any masked channel stops the extension
+   void searchForWings(std::list<std::pair<int, int> > &newlines,
+                           const casa::Vector<casa::Int> &signs)
+                           throw(casa::AipsError);
                            
    // An auxiliary object function to test whether two lines have a non-void