source: trunk/python/asaplinefind.py@ 420

Last change on this file since 420 was 371, checked in by vor010, 20 years ago

LineFinder/baseline fitter full functionality
auto_poly_baseline can automatically identify lines and fit baselines afterwards

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.7 KB
Line 
1import _asap
2
3class linefinder:
4 """
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
33 """
34
35 def __init__(self):
36 """
37 Create a line finder object.
38 """
39 self.finder = _asap.linefinder()
40 return
41
42 def set_options(self,threshold=1.7320508075688772,min_nchan=3,
43 avg_limit=8,box_size=0.2):
44 """
45 Set the parameters of the algorithm
46 Parameters:
47 threshold a single channel S/N ratio above which the
48 channel is considered to be a detection
49 Default is sqrt(3), which together with
50 min_nchan=3 gives a 3-sigma criterion
51 min_nchan a minimal number of consequtive channels,
52 which should satisfy a threshold criterion to
53 be a detection. Default is 3.
54 avg_limit A number of consequtive channels not greater than
55 this parameter can be averaged to search for
56 broad lines. Default is 8.
57 box_size A running mean box size specified as a fraction
58 of the total spectrum length. Default is 1/5
59 Note: For bad baselines threshold should be increased,
60 and avg_limit decreased (or even switched off completely by
61 setting this parameter to 1) to avoid detecting baseline
62 undulations instead of real lines.
63 """
64 self.finder.setoptions(threshold,min_nchan,avg_limit,box_size)
65 return
66
67 def set_scan(self,scan,mask=None,edge=(0,0)):
68 """
69 Set the 'data' (scantable) to work with.
70 Parameters:
71 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
75 specified, the same number will be dropped from
76 both sides of the spectrum. Default is to keep
77 all channels
78 """
79 if not scan:
80 raise RuntimeError, 'Please give a correct scan'
81 if len(edge)>2:
82 raise RuntimeError, "The edge parameter should have two \
83 or less elements"
84 if mask is None:
85 from numarray import ones
86 self.finder.setscan(scan,ones(scan.nchan()),edge)
87 else:
88 self.finder.setscan(scan,mask,edge)
89 return
90 def find_lines(self,nRow=0):
91 """
92 Search for spectral lines in the scan assigned in set_scan.
93 Current Beam/IF/Pol is used, Row is specified by parameter
94 A number of lines found will be returned
95 """
96 return self.finder.findlines(nRow)
97 def get_mask(self,invert=False):
98 """
99 Get the mask to mask out all lines that have been found (default)
100
101 Parameters:
102 invert if True, only channels belong to lines will be unmasked
103
104 Note: all channels originally masked by the input mask or
105 dropped out by the edge parameter will still be excluded
106 regardless on the invert option
107 """
108 return self.finder.getmask(invert)
109 def get_ranges(self,defunits=True):
110 """
111 Get ranges (start and end channels or velocities) for all spectral
112 lines found.
113
114 Parameters:
115 defunits if True (default), the range will use the same units
116 as set for the scan (e.g. LSR velocity)
117 if False, the range will be expressed in channels
118 """
119 if (defunits):
120 return self.finder.getlineranges()
121 else:
122 return self.finder.getlinerangesinchannels()
123
124def 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 workscan.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 x=f.get_parameters()
192 workscan._setspectrum(f.fitter.getresidual(),iRow)
193 workscan.set_cursor(sel[0],sel[1],sel[2])
194 workscan._vb = vb
195 if not insitu:
196 return workscan
Note: See TracBrowser for help on using the repository browser.