Changeset 370
- Timestamp:
- 02/06/05 19:46:16 (20 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/asaplinefind.py
r369 r370 4 4 """ 5 5 The class for automated spectral line search in ASAP. 6 7 Example: 8 fl=linefinder() 9 fl.set_scan(sc,edge=(50,)) 10 fl.set_options(threshold=3) 11 nlines=fl.find_lines() 12 if nlines!=0: 13 print "Found ",nlines," spectral lines" 14 print fl.get_ranges(False) 15 else: 16 print "No lines found!" 17 sc2=poly_baseline(sc,fl.get_mask(),7) 18 19 The algorithm involves a simple threshold criterion. The line is 20 considered to be detected if a specified number of consequtive 21 channels (default is 3) is brighter (with respect to the current baseline 22 estimate) than the threshold times the noise level. This criterion is 23 applied in the iterative procedure updating baseline estimate and trying 24 reduced spectral resolutions to detect broad lines as well. The off-line 25 noise level is determined at each iteration as an average of 80% of the 26 lowest variances across the spectrum (i.e. histogram equalization is 27 used to avoid missing weak lines if strong ones are present). For 28 bad baseline shapes it is reccommended to increase the threshold and 29 possibly switch the averaging option off (see set_options) to 30 detect strong lines only, fit a high order baseline and repeat the line 31 search. 32 6 33 """ 7 34 … … 61 88 self.finder.setscan(scan,mask,edge) 62 89 return 63 def find_lines(self ):90 def find_lines(self,nRow=0): 64 91 """ 65 92 Search for spectral lines in the scan assigned in set_scan. 93 Current Beam/IF/Pol is used, Row is specified by parameter 66 94 A number of lines found will be returned 67 95 """ 68 return self.finder.findlines( )96 return self.finder.findlines(nRow) 69 97 def get_mask(self,invert=False): 70 98 """ … … 89 117 if False, the range will be expressed in channels 90 118 """ 91 return self.finder.getlineranges(defunits) 119 if (defunits): 120 return self.finder.getlineranges() 121 else: 122 return self.finder.getlinerangesinchannels() 123 124 def auto_poly_baseline(scan, mask=None, edge=(0,0), order=0, 125 threshold=3,insitu=None): 126 """ 127 Return a scan which has been baselined (all rows) by a polynomial. 128 Spectral lines are detected first using linefinder and masked out 129 to avoid them affecting the baseline solution. 130 131 Parameters: 132 scan: a scantable 133 mask: an optional mask retreived from scantable 134 edge: an optional number of channel to drop at 135 the edge of spectrum. If only one value is 136 specified, the same number will be dropped from 137 both sides of the spectrum. Default is to keep 138 all channels 139 order: the order of the polynomial (default is 0) 140 threshold: the threshold used by line finder. It is better to 141 keep it large as only strong lines affect the 142 baseline solution. 143 insitu: if False a new scantable is returned. 144 Otherwise, the scaling is done in-situ 145 The default is taken from .asaprc (False) 146 147 Example: 148 sc2=auto_poly_baseline(sc,order=7) 149 """ 150 from asap.asapfitter import fitter 151 from asap import scantable 152 153 # setup fitter 154 155 f = fitter() 156 f._verbose(True) 157 f.set_function(poly=order) 158 159 # setup line finder 160 161 fl=linefinder() 162 fl.set_options(threshold=threshold) 163 164 if not insitu: 165 workscan=scan.copy() 166 else: 167 workscan=scan 168 169 vb=workscan._vb 170 # remember the verbose parameter and selection 171 workscan._vb=False 172 sel=workscan.get_cursor() 173 rows=range(workscan.nrow()) 174 for i in range(workscan.nbeam()): 175 workscan.setbeam(i) 176 for j in range(workscan.nif()): 177 workscan.setif(j) 178 for k in range(workscan.npol()): 179 scan.setpol(k) 180 if f._vb: 181 print "Processing:" 182 print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k) 183 for iRow in rows: 184 fl.set_scan(workscan,mask,edge) 185 fl.find_lines(iRow) 186 f.set_scan(workscan, fl.get_mask()) 187 f.x=workscan._getabcissa(iRow) 188 f.y=workscan._getspectrum(iRow) 189 f.data=None 190 f.fit() 191 workscan._setspectrum(f.getresidual(),iRow) 192 workscan.set_cursor(sel[0],sel[1],sel[2]) 193 workscan._vb = vb 194 if not insitu: 195 return scan -
trunk/src/SDLineFinder.cc
r369 r370 673 673 674 674 // search for spectral lines. Number of lines found is returned 675 int SDLineFinder::findLines( ) throw(casa::AipsError)675 int SDLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError) 676 676 { 677 677 const int minboxnchan=4; … … 685 685 throw AipsError("SDLineFinder::findLines - box_size is too small"); 686 686 687 scan->getSpectrum(spectrum );687 scan->getSpectrum(spectrum, whichRow); 688 688 689 689 lines.resize(0); // search from the scratch 690 last_row_used=whichRow; 690 691 Vector<Bool> temp_mask(mask); 691 692 … … 828 829 } 829 830 830 // get range for all lines found. If defunits is true (default), the 831 // same units as used in the scan will be returned (e.g. velocity 832 // instead of channels). If defunits is false, channels will be returned 833 std::vector<int> SDLineFinder::getLineRanges(bool defunits) 831 // get range for all lines found. The same units as used in the scan 832 // will be returned (e.g. velocity instead of channels). 833 std::vector<double> SDLineFinder::getLineRanges() 834 834 const throw(casa::AipsError) 835 { 836 // convert to required abscissa units 837 std::vector<double> vel=scan->getAbcissa(last_row_used); 838 std::vector<int> ranges=getLineRangesInChannels(); 839 std::vector<double> res(ranges.size()); 840 841 std::vector<int>::const_iterator cri=ranges.begin(); 842 std::vector<double>::iterator outi=res.begin(); 843 for (;cri!=ranges.end() && outi!=res.end();++cri,++outi) 844 if (uInt(*cri)>=vel.size()) 845 throw AipsError("SDLineFinder::getLineRanges - getAbcissa provided less channels than reqired"); 846 else *outi=vel[*cri]; 847 return res; 848 } 849 850 // The same as getLineRanges, but channels are always used to specify 851 // the range 852 std::vector<int> SDLineFinder::getLineRangesInChannels() 853 const throw(casa::AipsError) 835 854 { 836 855 try { 837 856 if (scan.null()) 838 throw AipsError("SDLineFinder::getLineRanges - a scan should be set first,"857 throw AipsError("SDLineFinder::getLineRangesInChannels - a scan should be set first," 839 858 " use set_scan followed by find_lines"); 840 859 DebugAssert(mask.nelements()==scan->nChan(), AipsError); 841 860 842 861 if (!lines.size()) 843 throw AipsError("SDLineFinder::getLineRanges - one have to search for "862 throw AipsError("SDLineFinder::getLineRangesInChannels - one have to search for " 844 863 "lines first, use find_lines"); 845 864 846 // temporary847 if (defunits)848 throw AipsError("SDLineFinder::getLineRanges - sorry, defunits=true have not "849 "yet been implemented");850 //851 865 std::vector<int> res(2*lines.size()); 852 866 // iterator through lines & result … … 857 871 if (++ri!=res.end()) 858 872 *ri=cli->second-1; 859 } 873 } 860 874 return res; 861 875 } … … 867 881 } 868 882 } 883 884 869 885 870 886 // an auxiliary function to remove all lines from the list, except the -
trunk/src/SDLineFinder.h
r369 r370 169 169 const boost::python::tuple &in_edge) throw(casa::AipsError); 170 170 171 // search for spectral lines. Number of lines found is returned 172 int findLines() throw(casa::AipsError); 171 // search for spectral lines for a row specified by whichRow and 172 // Beam/IF/Pol specified by current cursor set for the scantable 173 // Number of lines found is returned 174 int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError); 173 175 174 176 // 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); 180 182 181 // get range for all lines found. If defunits is true (default), the 182 // same units as used in the scan will be returned (e.g. velocity 183 // instead of channels). If defunits is false, channels will be returned 184 std::vector<int> getLineRanges(bool defunits=true) 185 const throw(casa::AipsError); 183 // get range for all lines found. The same units as used in the scan 184 // will be returned (e.g. velocity instead of channels). 185 std::vector<double> getLineRanges() const throw(casa::AipsError); 186 // The same as getLineRanges, but channels are always used to specify 187 // the range 188 std::vector<int> getLineRangesInChannels() const throw(casa::AipsError); 186 189 protected: 187 190 // auxiliary function to average adjacent channels and update the mask … … 236 239 // adjacent channels to search 237 240 // for broad lines. see setOptions 241 casa::uInt last_row_used; // the Row number specified 242 // during the last findLines call 238 243 std::list<std::pair<int, int> > lines; // container of start and stop+1 239 244 // channels of the spectral lines -
trunk/src/python_SDLineFinder.cc
r369 r370 45 45 .def("getmask",&SDLineFinder::getMask) 46 46 .def("getlineranges",&SDLineFinder::getLineRanges) 47 .def("getlinerangesinchannels",&SDLineFinder::getLineRangesInChannels) 47 48 ; 48 49 };
Note:
See TracChangeset
for help on using the changeset viewer.