Changeset 907


Ignore:
Timestamp:
03/20/06 15:17:06 (18 years ago)
Author:
vor010
Message:

LineFinder? & auto_poly_baseline: a support of
new scantable format has been added. Now findLines accept a mask and edge
parameters, which can therefore be different for different rows. Constructor
need now just a scan table. Boost types removed from STLineFinder. auto_poly_baseline can accept a nested tuple of edges, which would be interpreted as
different edge parameters for different IFs

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asaplinefind.py

    r889 r907  
    77    Example:
    88       fl=linefinder()
    9        fl.set_scan(sc,edge=(50,))
     9       fl.set_scan(sc)
    1010       fl.set_options(threshold=3)
    11        nlines=fl.find_lines()
     11       nlines=fl.find_lines(edge=(50,0))
    1212       if nlines!=0:
    1313          print "Found ",nlines," spectral lines"
    … …  
    6565        return
    6666
    67     def set_scan(self, scan, mask=None, edge=(0,0)):
     67    def set_scan(self, scan):
    6868        """
    6969        Set the 'data' (scantable) to work with.
    7070        Parameters:
    7171             scan:    a scantable
    72              mask:       an optional mask retreived from scantable
    73              edge:       an optional number of channel to drop at
    74                          the edge of spectrum. If only one value is
     72        """
     73        if not scan:
     74           raise RuntimeError, 'Please give a correct scan'
     75        return
     76    def find_lines(self,nRow=0,mask=None,edge=(0,0)):
     77        """
     78        Search for spectral lines in the scan assigned in set_scan.
     79        Parameters:
     80             nRow:       a row in the scantable to work with
     81             mask:       an optional mask (e.g. retreived from scantable)
     82             edge:       an optional number of channels to drop at
     83                         the edge of the spectrum. If only one value is
    7584                         specified, the same number will be dropped from
    7685                         both sides of the spectrum. Default is to keep
    7786                         all channels
     87        A number of lines found will be returned
    7888        """
    79         if not scan:
    80            raise RuntimeError, 'Please give a correct scan'
    81         if not scan._check_ifs():
    82            raise RuntimeError, 'IFs with different numbers of channels are not yet supported'
    83 
    8489        if isinstance(edge,int):
    8590           edge=(edge,)
    … …  
    96101        if mask is None:
    97102            from numarray import ones
    98             self.finder.setscan(scan,ones(scan.nchan(-1)),tuple(edge))
     103            return self.finder.findlines(ones(scan.nchan(nRow)),list(edge),nRow)
    99104        else:
    100             self.finder.setscan(scan,mask,tuple(edge))
    101         return
    102     def find_lines(self,nRow=0):
    103         """
    104         Search for spectral lines in the scan assigned in set_scan.
    105         Current Beam/IF/Pol is used, Row is specified by parameter
    106         A number of lines found will be returned
    107         """
    108         return self.finder.findlines(nRow)
     105            return self.finder.setscan(mask,list(edge),nRow)
    109106    def get_mask(self,invert=False):
    110107        """

  • trunk/python/scantable.py

    r889 r907  
    10331033                        specified, the same number will be dropped from
    10341034                        both sides of the spectrum. Default is to keep
    1035                         all channels
     1035                        all channels. Nested tuples represent individual
     1036                        edge selection for different IFs (a number of spectral
     1037                        channels can be different)
    10361038            order:      the order of the polynomial (default is 0)
    10371039            threshold:  the threshold used by line finder. It is better to
    … …  
    10511053        from asap import _is_sequence_or_number as _is_valid
    10521054
    1053         if not _is_valid(edge, int):
     1055        # check whether edge is set up for each IF individually
     1056        individualEdge = False;
     1057        if len(edge)>1:
     1058           if isinstance(edge[0],list) or isinstance(edge[0],tuple):
     1059           individualEdge = True;
     1060
     1061        if not _is_valid(edge, int) and not individualEdge:
    10541062            raise RuntimeError, "Parameter 'edge' has to be an integer or a \
    1055             pair of integers specified as a tuple"
     1063            pair of integers specified as a tuple. Nested tuples are allowed \
     1064            to make individual selection for different IFs."
     1065       
     1066        curedge = (0,0)
     1067        if individualEdge:
     1068           for edge_par in edge:
     1069               if not _is_valid(edge,int):
     1070                  raise RuntimeError, "Each element of the 'edge' tuple has \
     1071                  to be a pair of integers or an integer."
     1072        else:
     1073           curedge = edge;
    10561074
    10571075        # setup fitter
    … …  
    10671085        else:
    10681086            workscan=self
     1087
     1088        fl.set_scan(workscan)
    10691089
    10701090        rows=range(workscan.nrow())
    … …  
    10741094            msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" %        (workscan.getscan(r),workscan.getbeam(r),workscan.getif(r),workscan.getpol(r), workscan.getcycle(r))
    10751095            asaplog.push(msg, False)
    1076             fl.set_scan(workscan, mask, edge)
    1077             fl.find_lines(r)
     1096
     1097            # figure out edge parameter
     1098            if individualEdge:
     1099               if len(edge)>=workscan.getif(r):
     1100                  raise RuntimeError, "Number of edge elements appear to be less than the number of IFs"
     1101                  curedge = edge[workscan.getif(r)]
     1102           
     1103            # setup line finder
     1104            fl.find_lines(r,mask,curedge)
    10781105            f.set_scan(workscan, fl.get_mask())
    10791106            f.x = workscan._getabcissa(r)

  • trunk/src/Makefile

    r906 r907  
    1 TARGET   := /tmp/asap/_asap.so
     1TARGET   := ../tmp/asap/_asap.so
    22
    33# the casa environment AIPSPATH has to be defined
    … …  
    6464# python only 2.3 has been tested
    6565PYVERSION := 2.3
    66 PYTHONROOT := /usr
     66#PYTHONROOT := /usr
     67PYTHONROOT := /export/DEBIANlocal
    6768PYTHONINC := -I$(PYTHONROOT)/include/python$(PYVERSION)
    68 PYTHONLIB := -L$(PYTHONROOT)/lib -lpython$(PYVERSION)
     69PYTHONLIB := -L$(PYTHONROOT)/lib -lpython$(PYVERSION) -L/export/DEBIANlocal/gnu/lib/
    6970
    7071# has to be build with same g++ version as casa

  • trunk/src/STLineFinder.cpp

    r894 r907  
    150150                                           // the value - current mean
    151151                                           // (used to search wings)
     152   casa::Int last_sign;                    // a sign (+1, -1 or 0) of the
     153                                           // last point of the detected line
     154                                           //
    152155public:
    153156
    … …  
    189192   void processCurLine(const casa::Vector<casa::Bool> &mask)
    190193                                                 throw(casa::AipsError);
     194   
     195   // get the sign of runningBox->aboveMean(). The RunningBox pointer
     196   // should be defined
     197   casa::Int getAboveMeanSign() const throw();
    191198};
    192199
    … …  
    370377}
    371378
     379// get the sign of runningBox->aboveMean(). The RunningBox pointer
     380// should be defined
     381casa::Int LFAboveThreshold::getAboveMeanSign() const throw()
     382{
     383  const Float buf=running_box->aboveMean();
     384  if (buf>0) return 1;
     385  if (buf<0) return -1;
     386  return 0;
     387}
     388
    372389
    373390// process a channel: update cur_line and is_detected before and
    … …  
    377394{
    378395  try {
    379        if (detect) {
    380            if (is_detected_before)
    381                cur_line.second=running_box->getChannel()+1;
    382            else {
    383                is_detected_before=True;
    384                cur_line.first=running_box->getChannel();
    385                cur_line.second=running_box->getChannel()+1;
    386            }
    387        } else processCurLine(mask);
     396       if (is_detected_before) {
     397          // we have to check that the current detection has the
     398          // same sign of running_box->aboveMean
     399          // otherwise it could be a spurious detection
     400          if (last_sign && last_sign!=getAboveMeanSign())
     401              detect=False;
     402          }
     403          if (detect) {
     404              last_sign=getAboveMeanSign();
     405              if (is_detected_before)
     406                  cur_line.second=running_box->getChannel()+1;
     407              else {
     408                  is_detected_before=True;
     409                  cur_line.first=running_box->getChannel();
     410                  cur_line.second=running_box->getChannel()+1;
     411              }
     412          } else processCurLine(mask);
    388413  }
    389414  catch (const AipsError &ae) {
    … …  
    490515                  threshold*offline_variance), mask);
    491516           else processCurLine(mask); // just finish what was accumulated before
     517
     518           signs[ch]=getAboveMeanSign();
     519           // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
     520           //   threshold*offline_variance<<endl;
     521           
    492522           const Float buf=running_box->aboveMean();
    493523           if (buf>0) signs[ch]=1;
    … …  
    635665STLineFinder::~STLineFinder() throw(AipsError) {}
    636666
    637 // set scan to work with (in_scan parameter), associated mask (in_mask
    638 // parameter) and the edge channel rejection (in_edge parameter)
     667// set scan to work with (in_scan parameter)
     668void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
     669{
     670  scan=in_scan.getCP();
     671  AlwaysAssert(!scan.null(),AipsError);
     672 
     673}
     674
     675// search for spectral lines. Number of lines found is returned
     676// in_edge and in_mask control channel rejection for a given row
    639677//   if in_edge has zero length, all channels chosen by mask will be used
    640678//   if in_edge has one element only, it represents the number of
    … …  
    642680//   in_edge is introduced for convinience, although all functionality
    643681//   can be achieved using a spectrum mask only
    644 void STLineFinder::setScan(const ScantableWrapper &in_scan,
    645                const std::vector<bool> &in_mask,
    646                const boost::python::tuple &in_edge) throw(AipsError)
    647 {
    648   try {
    649        scan=in_scan.getCP();
    650        AlwaysAssert(!scan.null(),AipsError);
    651 
    652        mask=in_mask;
    653        if (mask.nelements()!=scan->nchan())
    654            throw AipsError("STLineFinder::setScan - in_scan and in_mask have different"
    655                            "number of spectral channels.");
    656 
    657        // number of elements in the in_edge tuple
    658        int n=extract<int>(in_edge.attr("__len__")());
    659        if (n>2 || n<0)
    660            throw AipsError("STLineFinder::setScan - the length of the in_edge parameter"
    661                            "should not exceed 2");
    662        if (!n) {
    663            // all spectra, no rejection
    664            edge.first=0;
    665            edge.second=scan->nchan();
    666        } else {
    667            edge.first=extract<int>(in_edge.attr("__getitem__")(0));
    668            if (edge.first<0)
    669                throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
    670                                 "number of channels to drop");
    671            if (edge.first>=scan->nchan())
    672                throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
    673            if (n==2) {
    674                edge.second=extract<int>(in_edge.attr("__getitem__")(1));
    675                if (edge.second<0)
    676                    throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
    677                                    "number of channels to drop");
    678                edge.second=scan->nchan()-edge.second;
    679            } else edge.second=scan->nchan()-edge.first;
    680            if (edge.second<0 || (edge.first>=edge.second))
    681                throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
    682        }
    683   }
    684   catch(const AipsError &ae) {
    685        // setScan is unsuccessfull, reset scan/mask/edge
    686        scan=CountedConstPtr<Scantable>(); // null pointer
    687        mask.resize(0);
    688        edge=pair<int,int>(0,0);
    689        throw;
    690   }
    691 }
    692 
    693 // search for spectral lines. Number of lines found is returned
    694 int STLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
     682int STLineFinder::findLines(const std::vector<bool> &in_mask,
     683                const std::vector<int> &in_edge,
     684                const casa::uInt &whichRow) throw(casa::AipsError)
    695685{
    696686  const int minboxnchan=4;
    … …  
    698688      throw AipsError("STLineFinder::findLines - a scan should be set first,"
    699689                      " use set_scan");
    700   DebugAssert(mask.nelements()==scan->nchan(), AipsError);
    701   int max_box_nchan=int(scan->nchan()*box_size); // number of channels in running
     690 
     691  // set up mask and edge rejection
     692  mask=in_mask;
     693  if (mask.nelements()!=scan->nchan(whichRow))
     694      throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
     695            "number of spectral channels.");
     696
     697  // number of elements in in_edge
     698  if (in_edge.size()>2)
     699      throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
     700                      "should not exceed 2");
     701      if (!in_edge.size()) {
     702           // all spectra, no rejection
     703           edge.first=0;
     704           edge.second=scan->nchan(whichRow);
     705      } else {
     706           edge.first=in_edge[0];
     707           if (edge.first<0)
     708               throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
     709                                "number of channels to drop");
     710           if (edge.first>=scan->nchan(whichRow))
     711               throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
     712           if (in_edge.size()==2) {
     713               edge.second=in_edge[1];
     714               if (edge.second<0)
     715                   throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
     716                                   "number of channels to drop");
     717               edge.second=scan->nchan(whichRow)-edge.second;
     718           } else edge.second=scan->nchan(whichRow)-edge.first;
     719           if (edge.second<0 || (edge.first>=edge.second))
     720               throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
     721       }
     722     
     723  //
     724  int max_box_nchan=int(scan->nchan(whichRow)*box_size); // number of channels in running
    702725                                                 // box
    703726  if (max_box_nchan<2)
    … …  
    833856           throw AipsError("STLineFinder::getMask - a scan should be set first,"
    834857                      " use set_scan followed by find_lines");
    835        DebugAssert(mask.nelements()==scan->nchan(), AipsError);
     858       DebugAssert(mask.nelements()==scan->nchan(last_row_used), AipsError);
    836859       /*
    837860       if (!lines.size())
    … …  
    892915           throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
    893916                      " use set_scan followed by find_lines");
    894        DebugAssert(mask.nelements()==scan->nchan(), AipsError);
     917       DebugAssert(mask.nelements()==scan->nchan(last_row_used), AipsError);
    895918
    896919       if (!lines.size())

  • trunk/src/STLineFinder.h

    r891 r907  
    158158                   const casa::Float &in_box_size=0.2) throw();
    159159
    160    // set the scan to work with (in_scan parameter), associated mask (in_mask
    161    // parameter) and the edge channel rejection (in_edge parameter)
    162    //   if in_edge has zero length, all channels chosen by mask will be used
     160   // set the scan to work with (in_scan parameter)
     161   void setScan(const ScantableWrapper &in_scan) throw(casa::AipsError);
     162
     163   // search for spectral lines in a row specified by whichRow
     164   // in_mask and in_edge parameters control channel rejection
     165   //  if in_edge has zero length, all channels chosen by mask will be used
    163166   //   if in_edge has one element only, it represents the number of
    164167   //      channels to drop from both sides of the spectrum
    165168   //   in_edge is introduced for convinience, although all functionality
    166169   //   can be achieved using a spectrum mask only
    167    void setScan(const ScantableWrapper &in_scan,
    168                 const std::vector<bool> &in_mask,
    169                 const boost::python::tuple &in_edge) throw(casa::AipsError);
    170 
    171    // search for spectral lines for a row specified by whichRow and
    172    // Beam/IF/Pol specified by current cursor set for the scantable
    173170   // Number of lines found is returned
    174    int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError);
     171   int findLines(const std::vector<bool> &in_mask,
     172                 const std::vector<int> &in_edge = std::vector<int>(),
     173                 const casa::uInt &whichRow = 0) throw(casa::AipsError);
    175174
    176175   // get the mask to mask out all lines that have been found (default)
Note: See TracChangeset for help on using the changeset viewer.