Changeset 370 for trunk

Timestamp:

02/06/05 19:46:16 (20 years ago)

Author:

vor010

Message:

LineFinder: help is improved. Initial code for
automatic baseline fitter is written (not debugged yet)

Location:

trunk

Files:

• python/asaplinefind.py (modified) (3 diffs)
• src/SDLineFinder.cc (modified) (5 diffs)
• src/SDLineFinder.h (modified) (3 diffs)
• src/python_SDLineFinder.cc (modified) (1 diff)

trunk/python/asaplinefind.py

@@ -4 +4 @@
     """
     The class for automated spectral line search in ASAP.
+
+    Example:
+       fl=linefinder()
+       fl.set_scan(sc,edge=(50,))
+       fl.set_options(threshold=3)
+       nlines=fl.find_lines()
+       if nlines!=0:
+          print "Found ",nlines," spectral lines"
+          print fl.get_ranges(False)
+       else:
+          print "No lines found!"
+       sc2=poly_baseline(sc,fl.get_mask(),7)
+    
+    The algorithm involves a simple threshold criterion. The line is
+    considered to be detected if a specified number of consequtive
+    channels (default is 3) is brighter (with respect to the current baseline
+    estimate) than the threshold times the noise level. This criterion is
+    applied in the iterative procedure updating baseline estimate and trying
+    reduced spectral resolutions to detect broad lines as well. The off-line
+    noise level is determined at each iteration as an average of 80% of the
+    lowest variances across the spectrum (i.e. histogram equalization is
+    used to avoid missing weak lines if strong ones are present). For
+    bad baseline shapes it is reccommended to increase the threshold and 
+    possibly switch the averaging option off (see set_options) to
+    detect strong lines only, fit a high order baseline and repeat the line
+    search. 
+
     """
 
@@ -61 +88 @@
             self.finder.setscan(scan,mask,edge)
         return 
-    def find_lines(self):
+    def find_lines(self,nRow=0):
         """
         Search for spectral lines in the scan assigned in set_scan.
+        Current Beam/IF/Pol is used, Row is specified by parameter
         A number of lines found will be returned
         """
-        return self.finder.findlines()
+        return self.finder.findlines(nRow)
     def get_mask(self,invert=False):
         """
@@ -89 +117 @@
                         if False, the range will be expressed in channels 
         """
-        return self.finder.getlineranges(defunits)
+        if (defunits):
+            return self.finder.getlineranges()
+        else:
+            return self.finder.getlinerangesinchannels()
+
+def auto_poly_baseline(scan, mask=None, edge=(0,0), order=0,
+    threshold=3,insitu=None):
+    """
+    Return a scan which has been baselined (all rows) by a polynomial.
+    Spectral lines are detected first using linefinder and masked out
+    to avoid them affecting the baseline solution.
+
+    Parameters:
+        scan:    a scantable
+        mask:       an optional mask retreived from scantable
+        edge:       an optional number of channel to drop at
+                    the edge of spectrum. If only one value is
+                    specified, the same number will be dropped from
+                    both sides of the spectrum. Default is to keep
+                    all channels
+        order:      the order of the polynomial (default is 0)
+        threshold:  the threshold used by line finder. It is better to
+                    keep it large as only strong lines affect the
+                    baseline solution.
+        insitu:     if False a new scantable is returned.
+                    Otherwise, the scaling is done in-situ
+                    The default is taken from .asaprc (False)
+
+    Example:
+        sc2=auto_poly_baseline(sc,order=7)
+    """
+    from asap.asapfitter import fitter
+    from asap import scantable
+
+    # setup fitter
+
+    f = fitter()
+    f._verbose(True)
+    f.set_function(poly=order)
+
+    # setup line finder
+
+    fl=linefinder()
+    fl.set_options(threshold=threshold)
+
+    if not insitu:
+        workscan=scan.copy()
+    else:
+        workscan=scan
+
+    vb=workscan._vb
+    # remember the verbose parameter and selection
+    workscan._vb=False
+    sel=workscan.get_cursor()
+    rows=range(workscan.nrow()) 
+    for i in range(workscan.nbeam()):
+        workscan.setbeam(i)
+        for j in range(workscan.nif()):
+            workscan.setif(j)
+            for k in range(workscan.npol()):
+                scan.setpol(k)
+                if f._vb:
+                   print "Processing:"
+                   print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
+                for iRow in rows:
+                   fl.set_scan(workscan,mask,edge)
+                   fl.find_lines(iRow)
+                   f.set_scan(workscan, fl.get_mask())
+                   f.x=workscan._getabcissa(iRow)
+                   f.y=workscan._getspectrum(iRow)
+                   f.data=None
+                   f.fit()
+                   workscan._setspectrum(f.getresidual(),iRow)
+    workscan.set_cursor(sel[0],sel[1],sel[2])
+    workscan._vb = vb
+    if not insitu:
+       return scan

trunk/src/SDLineFinder.cc

@@ -673 +673 @@
 
 // search for spectral lines. Number of lines found is returned
-int SDLineFinder::findLines() throw(casa::AipsError)
+int SDLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
 { 
   const int minboxnchan=4;
@@ -685 +685 @@
       throw AipsError("SDLineFinder::findLines - box_size is too small");
 
-  scan->getSpectrum(spectrum);
+  scan->getSpectrum(spectrum, whichRow);
 
   lines.resize(0); // search from the scratch
+  last_row_used=whichRow;
   Vector<Bool> temp_mask(mask);
 
@@ -828 +829 @@
 }
 
-// get range for all lines found. If defunits is true (default), the
-// same units as used in the scan will be returned (e.g. velocity
-// instead of channels). If defunits is false, channels will be returned
-std::vector<int> SDLineFinder::getLineRanges(bool defunits)
+// get range for all lines found. The same units as used in the scan
+// will be returned (e.g. velocity instead of channels).
+std::vector<double> SDLineFinder::getLineRanges()
                              const throw(casa::AipsError)
+{
+  // convert to required abscissa units
+  std::vector<double> vel=scan->getAbcissa(last_row_used);
+  std::vector<int> ranges=getLineRangesInChannels();
+  std::vector<double> res(ranges.size());
+
+  std::vector<int>::const_iterator cri=ranges.begin();
+  std::vector<double>::iterator outi=res.begin();
+  for (;cri!=ranges.end() && outi!=res.end();++cri,++outi)
+       if (uInt(*cri)>=vel.size())
+          throw AipsError("SDLineFinder::getLineRanges - getAbcissa provided less channels than reqired");
+       else *outi=vel[*cri];
+  return res;
+}
+
+// The same as getLineRanges, but channels are always used to specify
+// the range
+std::vector<int> SDLineFinder::getLineRangesInChannels()
+                                   const throw(casa::AipsError)
 {
   try {
        if (scan.null())
-           throw AipsError("SDLineFinder::getLineRanges - a scan should be set first,"
+           throw AipsError("SDLineFinder::getLineRangesInChannels - a scan should be set first,"
                       " use set_scan followed by find_lines");
        DebugAssert(mask.nelements()==scan->nChan(), AipsError);
        
        if (!lines.size())
-           throw AipsError("SDLineFinder::getLineRanges - one have to search for "
+           throw AipsError("SDLineFinder::getLineRangesInChannels - one have to search for "
                            "lines first, use find_lines");
                            
-       // temporary
-       if (defunits)
-           throw AipsError("SDLineFinder::getLineRanges - sorry, defunits=true have not "
-                           "yet been implemented");
-       //
        std::vector<int> res(2*lines.size());
        // iterator through lines & result
@@ -857 +871 @@
             if (++ri!=res.end()) 
                 *ri=cli->second-1;          
-       }
+       }       
        return res;
   }
@@ -867 +881 @@
   }
 }
+
+
 
 // an auxiliary function to remove all lines from the list, except the

trunk/src/SDLineFinder.h

@@ -169 +169 @@
                 const boost::python::tuple &in_edge) throw(casa::AipsError);
 
-   // search for spectral lines. Number of lines found is returned
-   int findLines() throw(casa::AipsError);
+   // search for spectral lines for a row specified by whichRow and
+   // Beam/IF/Pol specified by current cursor set for the scantable
+   // Number of lines found is returned   
+   int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError);
 
    // get the mask to mask out all lines that have been found (default)
@@ -179 +181 @@
    std::vector<bool> getMask(bool invert=false) const throw(casa::AipsError);
 
-   // get range for all lines found. If defunits is true (default), the
-   // same units as used in the scan will be returned (e.g. velocity
-   // instead of channels). If defunits is false, channels will be returned
-   std::vector<int>   getLineRanges(bool defunits=true)
-                                const throw(casa::AipsError);
+   // get range for all lines found. The same units as used in the scan
+   // will be returned (e.g. velocity instead of channels).   
+   std::vector<double>   getLineRanges() const throw(casa::AipsError);
+   // The same as getLineRanges, but channels are always used to specify
+   // the range
+   std::vector<int> getLineRangesInChannels() const throw(casa::AipsError);
 protected:
    // auxiliary function to average adjacent channels and update the mask
@@ -236 +239 @@
                                            // adjacent channels to search
                                            // for broad lines. see setOptions
+   casa::uInt last_row_used;                // the Row number specified
+                                           // during the last findLines call
    std::list<std::pair<int, int> > lines;  // container of start and stop+1
                                            // channels of the spectral lines

trunk/src/python_SDLineFinder.cc

@@ -45 +45 @@
          .def("getmask",&SDLineFinder::getMask)
          .def("getlineranges",&SDLineFinder::getLineRanges)
+         .def("getlinerangesinchannels",&SDLineFinder::getLineRangesInChannels)
        ;
      };