source: trunk/src/SDLineFinder.cc@ 323

Last change on this file since 323 was 297, checked in by vor010, 20 years ago

SDLineFinder: C++ class and python binder
have been added. This is an initial version, but it works in some simple

cases. Makefile and install.sh were updated to account for new source files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.7 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDLineFinder.cc: A class for automated spectral line search
3//#--------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32
33// ASAP
34#include "SDLineFinder.h"
35
36// STL
37#include <iostream>
38
39using namespace asap;
40using namespace casa;
41using namespace std;
42using namespace boost::python;
43
44// SDLineFinder - a class for automated spectral line search
45
46SDLineFinder::SDLineFinder() throw() : edge(0,0)
47{
48 // detection threshold - the minimal signal to noise ratio
49 threshold=3.; // 3 sigma is a default
50 box_size=1./16.; // default box size for running mean calculations is
51 // 1/16 of the whole spectrum
52 // A minimum number of consequtive channels, which should satisfy
53 // the detection criterion, to be a detection
54 min_nchan=3; // default is 3 channels
55}
56
57SDLineFinder::~SDLineFinder() throw(AipsError) {}
58
59// set scan to work with (in_scan parameter), associated mask (in_mask
60// parameter) and the edge channel rejection (in_edge parameter)
61// if in_edge has zero length, all channels chosen by mask will be used
62// if in_edge has one element only, it represents the number of
63// channels to drop from both sides of the spectrum
64// in_edge is introduced for convinience, although all functionality
65// can be achieved using a spectrum mask only
66void SDLineFinder::setScan(const SDMemTableWrapper &in_scan,
67 const std::vector<bool> &in_mask,
68 const boost::python::tuple &in_edge) throw(AipsError)
69{
70 try {
71 scan=in_scan.getCP();
72 AlwaysAssert(!scan.null(),AipsError);
73 if (scan->nRow()!=1)
74 throw AipsError("SDLineFinder::setScan - in_scan contains more than 1 row."
75 "Choose one first.");
76 mask=in_mask;
77 if (mask.nelements()!=scan->nChan())
78 throw AipsError("SDLineFinder::setScan - in_scan and in_mask have different"
79 "number of spectral channels.");
80
81 // number of elements in the in_edge tuple
82 int n=extract<int>(in_edge.attr("__len__")());
83 if (n>2 || n<0)
84 throw AipsError("SDLineFinder::setScan - the length of the in_edge parameter"
85 "should not exceed 2");
86 if (!n) {
87 // all spectrum, no rejection
88 edge.first=0;
89 edge.second=scan->nChan();
90 } else {
91 edge.first=extract<int>(in_edge.attr("__getitem__")(0));
92 if (edge.first<0)
93 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
94 "number of channels to drop");
95 if (edge.first>=scan->nChan())
96 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
97 if (n==2) {
98 edge.second=extract<int>(in_edge.attr("__getitem__")(1));
99 if (edge.second<0)
100 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
101 "number of channels to drop");
102 edge.second=scan->nChan()-edge.second;
103 } else edge.second=scan->nChan()-edge.first;
104 if (edge.second<0 || (edge.second+edge.first)>scan->nChan())
105 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
106 }
107 }
108 catch(const AipsError &ae) {
109 // setScan is unsuccessfull, reset scan/mask/edge
110 scan=CountedConstPtr<SDMemTable>(); // null pointer
111 mask.resize(0);
112 edge=pair<int,int>(0,0);
113 throw;
114 }
115}
116
117// search for spectral lines. Number of lines found is returned
118int SDLineFinder::findLines() throw(casa::AipsError)
119{
120 const int minboxnchan=4;
121 if (scan.null())
122 throw AipsError("SDLineFinder::findLines - a scan should be set first,"
123 " use set_scan");
124 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
125 lines.resize(0); // search from the scratch
126
127 max_box_nchan=int(scan->nChan()*box_size); // number of channels in running
128 // box
129
130 if (max_box_nchan<2)
131 throw AipsError("SDLineFinder::findLines - box_size is too small");
132
133 scan->getSpectrum(spectrum);
134
135 // fill statistics for initial box
136 box_chan_cntr=0; // no channels are currently in the box
137 sum=0; // initialize statistics
138 sumsq=0;
139 int initial_box_ch=edge.first;
140 for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan;
141 ++initial_box_ch)
142 advanceRunningBox(initial_box_ch);
143
144 if (initial_box_ch==edge.second)
145 throw AipsError("SDLineFinder::findLines - too much channels are masked");
146
147 // actual search algorithm
148 is_detected_before=False;
149
150 if (box_chan_cntr>=minboxnchan)
151 // there is a minimum amount of data. We can search in the
152 // half of the initial box
153 for (int n=edge.first;n<initial_box_ch-max_box_nchan/2;++n)
154 processChannel(n);
155
156 // now the box can be moved. n+max_box_nchan/2 is a new index which haven't
157 // yet been included in the running mean.
158 for (int n=initial_box_ch-max_box_nchan/2;n<edge.second;++n) {
159 advanceRunningBox(n+max_box_nchan/2); // update running mean & variance
160 if (box_chan_cntr>=minboxnchan) // have enough data to process
161 processChannel(n);
162 else processCurLine(); // just finish what was accumulated before
163 }
164 return int(lines.size());
165}
166
167// supplementary function to control running mean calculations.
168// It adds a specified channel to the running mean box and
169// removes (ch-max_box_nchan+1)'th channel from there
170// Channels, for which the mask is false or index is beyond the
171// allowed range, are ignored
172void SDLineFinder::advanceRunningBox(int ch) throw(AipsError)
173{
174 if (ch>=edge.first && ch<edge.second)
175 if (mask[ch]) { // ch is a valid channel
176 ++box_chan_cntr;
177 sum+=spectrum[ch];
178 sumsq+=square(spectrum[ch]);
179 }
180 int ch2remove=ch-max_box_nchan;
181 if (ch2remove>=edge.first && ch2remove<edge.second)
182 if (mask[ch2remove]) { // ch2remove is a valid channel
183 --box_chan_cntr;
184 sum-=spectrum[ch2remove];
185 sumsq-=square(spectrum[ch2remove]);
186 }
187}
188
189// test a channel against current running mean & rms
190// if channel specified is masked out or beyond the allowed indexes,
191// false is returned
192Bool SDLineFinder::testChannel(int ch) throw(exception, AipsError)
193{
194 if (ch<edge.first || ch>=edge.second) return False;
195 if (!mask[ch]) return False;
196 DebugAssert(box_chan_cntr, AipsError);
197 Float mean=sum/Float(box_chan_cntr);
198 Float variance=sqrt(sumsq/Float(box_chan_cntr)-square(mean));
199 /*
200 if (ch>3900 && ch<4100)
201 cout<<"Tested "<<ch<<" mean="<<mean<<" variance="<<variance<<" sp-mean="<<spectrum[ch]-mean<<endl;
202 */
203 return fabs(spectrum[ch]-mean)>=threshold*variance;
204}
205
206// process a channel: update cur_line and is_detected before and
207// add a new line to the list, if necessary
208void SDLineFinder::processChannel(int ch) throw(casa::AipsError)
209{
210 try {
211 if (testChannel(ch)) {
212 if (is_detected_before)
213 cur_line.second=ch+1;
214 else {
215 is_detected_before=True;
216 cur_line.first=ch;
217 cur_line.second=ch+1;
218 }
219 } else processCurLine();
220 }
221 catch (const AipsError &ae) {
222 throw;
223 }
224 catch (const exception &ex) {
225 throw AipsError(String("SDLineFinder::processChannel - STL error: ")+ex.what());
226 }
227}
228
229// process the interval of channels stored in cur_line
230// if it satisfies the criterion, add this interval as a new line
231void SDLineFinder::processCurLine() throw(casa::AipsError)
232{
233 try {
234 if (is_detected_before) {
235 if (cur_line.second-cur_line.first>min_nchan) {
236 // it was a detection. We need to change the list
237 Bool add_new_line=False;
238 if (lines.size()) {
239 for (int i=lines.back().second;i<cur_line.first;++i)
240 if (mask[i]) { // one valid channel in between
241 // means that we deal with a separate line
242 add_new_line=True;
243 break;
244 }
245 } else add_new_line=True;
246 if (add_new_line)
247 lines.push_back(cur_line);
248 else lines.back().second=cur_line.second;
249 }
250 is_detected_before=False;
251 }
252 }
253 catch (const AipsError &ae) {
254 throw;
255 }
256 catch (const exception &ex) {
257 throw AipsError(String("SDLineFinder::processCurLine - STL error: ")+ex.what());
258 }
259}
260
261// get the mask to mask out all lines that have been found (default)
262// if invert=true, only channels belong to lines will be unmasked
263// Note: all channels originally masked by the input mask (in_mask
264// in setScan) or dropped out by the edge parameter (in_edge
265// in setScan) are still excluded regardless on the invert option
266std::vector<bool> SDLineFinder::getMask(bool invert)
267 const throw(casa::AipsError)
268{
269 try {
270 if (scan.null())
271 throw AipsError("SDLineFinder::getMask - a scan should be set first,"
272 " use set_scan followed by find_lines");
273 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
274 /*
275 if (!lines.size())
276 throw AipsError("SDLineFinder::getMask - one have to search for "
277 "lines first, use find_lines");
278 */
279 std::vector<bool> res_mask(mask.nelements());
280 // iterator through lines
281 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
282 for (int ch=0;ch<res_mask.size();++ch)
283 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
284 else if (!mask[ch]) res_mask[ch]=false;
285 else {
286 res_mask[ch]=!invert; // no line by default
287 if (cli==lines.end()) continue;
288 if (ch>=cli->first && ch<cli->second)
289 res_mask[ch]=invert; // this is a line
290 if (ch>=cli->second)
291 ++cli; // next line in the list
292 }
293
294 return res_mask;
295 }
296 catch (const AipsError &ae) {
297 throw;
298 }
299 catch (const exception &ex) {
300 throw AipsError(String("SDLineFinder::getMask - STL error: ")+ex.what());
301 }
302}
303
304// get range for all lines found. If defunits is true (default), the
305// same units as used in the scan will be returned (e.g. velocity
306// instead of channels). If defunits is false, channels will be returned
307std::vector<int> SDLineFinder::getLineRanges(bool defunits)
308 const throw(casa::AipsError)
309{
310 try {
311 if (scan.null())
312 throw AipsError("SDLineFinder::getLineRanges - a scan should be set first,"
313 " use set_scan followed by find_lines");
314 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
315
316 if (!lines.size())
317 throw AipsError("SDLineFinder::getLineRanges - one have to search for "
318 "lines first, use find_lines");
319
320 // temporary
321 if (defunits)
322 throw AipsError("SDLineFinder::getLineRanges - sorry, defunits=true have not "
323 "yet been implemented");
324 //
325 std::vector<int> res(2*lines.size());
326 // iterator through lines & result
327 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
328 std::vector<int>::iterator ri=res.begin();
329 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
330 *ri=cli->first;
331 if (++ri!=res.end())
332 *ri=cli->second-1;
333 }
334 return res;
335 }
336 catch (const AipsError &ae) {
337 throw;
338 }
339 catch (const exception &ex) {
340 throw AipsError(String("SDLineFinder::getLineRanges - STL error: ")+ex.what());
341 }
342}
Note: See TracBrowser for help on using the repository browser.