source: branches/hpc34/src/STLineFinder.h@ 2956

Last change on this file since 2956 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# STLineFinder.h: 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: STLineFinder.h 2580 2012-06-28 04:22:10Z ShinnosukeKawakami $
30//#---------------------------------------------------------------------------
31#ifndef STLINEFINDER_H
32#define STLINEFINDER_H
33
34// STL
35#include <vector>
36#include <list>
37#include <utility>
38#include <exception>
39
40// AIPS++
41#include <casa/aips.h>
42#include <casa/Exceptions/Error.h>
43#include <casa/Arrays/Vector.h>
44#include <casa/Utilities/Assert.h>
45#include <casa/Utilities/CountedPtr.h>
46
47// ASAP
48#include "ScantableWrapper.h"
49#include "Scantable.h"
50#include "STFitter.h"
51
52namespace asap {
53
54///////////////////////////////////////////////////////////////////////////////
55//
56// LFLineListOperations - a class incapsulating operations with line lists
57// The LF prefix stands for Line Finder
58//
59
60struct LFLineListOperations {
61 // concatenate two lists preserving the order. If two lines appear to
62 // be adjacent or have a non-void intersection, they are joined into
63 // the new line
64 static void addNewSearchResult(const std::list<std::pair<int, int> >
65 &newlines, std::list<std::pair<int, int> > &lines_list)
66 throw(casa::AipsError);
67
68 // extend all line ranges to the point where a value stored in the
69 // specified vector changes (e.g. value-mean change its sign)
70 // This operation is necessary to include line wings, which are below
71 // the detection threshold. If lines becomes adjacent, they are
72 // merged together. Any masked channel stops the extension
73 static void searchForWings(std::list<std::pair<int, int> > &newlines,
74 const casa::Vector<casa::Int> &signs,
75 const casa::Vector<casa::Bool> &mask,
76 const std::pair<int,int> &edge)
77 throw(casa::AipsError);
78protected:
79
80 // An auxiliary object function to test whether two lines have a non-void
81 // intersection
82 class IntersectsWith : public std::unary_function<pair<int,int>, bool> {
83 std::pair<int,int> line1; // range of the first line
84 // start channel and stop+1
85 public:
86 explicit IntersectsWith(const std::pair<int,int> &in_line1);
87 // return true if line2 intersects with line1 with at least one
88 // common channel, and false otherwise
89 bool operator()(const std::pair<int,int> &line2) const throw();
90 };
91
92 // An auxiliary object function to build a union of several lines
93 // to account for a possibility of merging the nearby lines
94 class BuildUnion {
95 std::pair<int,int> temp_line; // range of the first line
96 // start channel and stop+1
97 public:
98 explicit BuildUnion(const std::pair<int,int> &line1);
99 // update temp_line with a union of temp_line and new_line
100 // provided there is no gap between the lines
101 void operator()(const std::pair<int,int> &new_line) throw();
102 // return the result (temp_line)
103 const std::pair<int,int>& result() const throw();
104 };
105
106 // An auxiliary object function to test whether a specified line
107 // is at lower spectral channels (to preserve the order in the line list)
108 class LaterThan : public std::unary_function<pair<int,int>, bool> {
109 std::pair<int,int> line1; // range of the first line
110 // start channel and stop+1
111 public:
112 explicit LaterThan(const std::pair<int,int> &in_line1);
113
114 // return true if line2 should be placed later than line1
115 // in the ordered list (so, it is at greater channel numbers)
116 bool operator()(const std::pair<int,int> &line2) const throw();
117 };
118
119
120};
121
122//
123///////////////////////////////////////////////////////////////////////////////
124
125///////////////////////////////////////////////////////////////////////////////
126//
127// STLineFinder - a class for automated spectral line search
128//
129//
130
131struct STLineFinder : protected LFLineListOperations {
132 STLineFinder() throw();
133 virtual ~STLineFinder() throw(casa::AipsError);
134
135 // set the parameters controlling algorithm
136 // in_threshold a single channel threshold default is sqrt(3), which
137 // means together with 3 minimum channels at least 3 sigma
138 // detection criterion
139 // For bad baseline shape, in_threshold may need to be
140 // increased
141 // in_min_nchan minimum number of channels above the threshold to report
142 // a detection, default is 3
143 // in_avg_limit perform the averaging of no more than in_avg_limit
144 // adjacent channels to search for broad lines
145 // Default is 8, but for a bad baseline shape this
146 // parameter should be decreased (may be even down to a
147 // minimum of 1 to disable this option) to avoid
148 // confusing of baseline undulations with a real line.
149 // Setting a very large value doesn't usually provide
150 // valid detections.
151 // in_box_size the box size for running mean calculation. Default is
152 // 1./5. of the whole spectrum size
153 // in_noise_box the box size for off-line noise estimation (if working with
154 // local noise. Negative value means use global noise estimate
155 // Default is -1 (i.e. estimate using the whole spectrum)
156 // in_median true if median statistics is used as opposed to average of
157 // the lowest 80% of deviations (default)
158 void setOptions(const casa::Float &in_threshold=sqrt(3.),
159 const casa::Int &in_min_nchan=3,
160 const casa::Int &in_avg_limit=8,
161 const casa::Float &in_box_size=0.2,
162 const casa::Float &in_noise_box=-1.,
163 const casa::Bool &in_median = casa::False) throw();
164
165 void setDetailedOptions( const casa::Int &order=9 ) ;
166
167 // set the scan to work with (in_scan parameter)
168 void setScan(const ScantableWrapper &in_scan) throw(casa::AipsError);
169
170 // set spectrum data to work with. this is a method to allow linefinder work
171 // without setting scantable for the purpose of using linefinder inside some
172 // method in scantable class. (Dec. 22, 2010 by W.Kawasaki)
173 void setData(const std::vector<float> &in_spectrum);
174
175 // search for spectral lines in a row specified by whichRow
176 // in_mask and in_edge parameters control channel rejection
177 // if in_edge has zero length, all channels chosen by mask will be used
178 // if in_edge has one element only, it represents the number of
179 // channels to drop from both sides of the spectrum
180 // in_edge is introduced for convinience, although all functionality
181 // can be achieved using a spectrum mask only
182 // Number of lines found is returned
183 int findLines(const std::vector<bool> &in_mask,
184 const std::vector<int> &in_edge = std::vector<int>(),
185 const casa::uInt &whichRow = 0) throw(casa::AipsError);
186
187 // get the mask to mask out all lines that have been found (default)
188 // if invert=true, only channels belong to lines will be unmasked
189 // Note: all channels originally masked by the input mask (in_mask
190 // in setScan) or dropped out by the edge parameter (in_edge
191 // in setScan) are still excluded regardless on the invert option
192 std::vector<bool> getMask(bool invert=false) const throw(casa::AipsError);
193
194 // get range for all lines found. The same units as used in the scan
195 // will be returned (e.g. velocity instead of channels).
196 std::vector<double> getLineRanges() const throw(casa::AipsError);
197 // The same as getLineRanges, but channels are always used to specify
198 // the range
199 std::vector<int> getLineRangesInChannels() const throw(casa::AipsError);
200protected:
201 // auxiliary function to average adjacent channels and update the mask
202 // if at least one channel involved in summation is masked, all
203 // output channels will be masked. This function works with the
204 // spectrum and edge fields of this class, but updates the mask
205 // array specified, rather than the field of this class
206 // boxsize - a number of adjacent channels to average
207 void averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update,
208 const casa::Int &boxsize)
209 throw(casa::AipsError);
210
211 // auxiliary function to fit and subtract a polynomial from the current
212 // spectrum. It uses the Fitter class. This action is required before
213 // reducing the spectral resolution if the baseline shape is bad
214 void subtractBaseline(const casa::Vector<casa::Bool> &temp_mask,
215 const casa::Int &order) throw(casa::AipsError);
216
217 // an auxiliary function to remove all lines from the list, except the
218 // strongest one (by absolute value). If the lines removed are real,
219 // they will be find again at the next iteration. This approach
220 // increases the number of iterations required, but is able to remove
221 // the sidelobes likely to occur near strong lines.
222 // Later a better criterion may be implemented, e.g.
223 // taking into consideration the brightness of different lines. Now
224 // use the simplest solution
225 // temp_mask - mask to work with (may be different from original mask as
226 // the lines previously found may be masked)
227 // lines2update - a list of lines to work with
228 // nothing will be done if it is empty
229 // max_box_nchan - channels in the running box for baseline filtering
230 void keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask,
231 std::list<std::pair<int, int> > &lines2update,
232 int max_box_nchan)
233 throw (casa::AipsError);
234private:
235 casa::CountedConstPtr<Scantable> scan; // the scan to work with
236 casa::Vector<casa::Bool> mask; // associated mask
237 std::pair<int,int> edge; // start and stop+1 channels
238 // to work with
239 casa::Float threshold; // detection threshold - the
240 // minimal signal to noise ratio
241 casa::Double box_size; // size of the box for running
242 // mean calculations, specified as
243 // a fraction of the whole spectrum
244 int min_nchan; // A minimum number of consequtive
245 // channels, which should satisfy
246 // the detection criterion, to be
247 // a detection
248 casa::Int avg_limit; // perform the averaging of no
249 // more than in_avg_limit
250 // adjacent channels to search
251 // for broad lines. see setOptions
252 casa::uInt last_row_used; // the Row number specified
253 // during the last findLines call
254 std::list<std::pair<int, int> > lines; // container of start and stop+1
255 // channels of the spectral lines
256 // a buffer for the spectrum
257 mutable casa::Vector<casa::Float> spectrum;
258
259 // the box size for off-line noise estimation (if working with
260 // local noise. Negative value means use global noise estimate
261 // Default is -1 (i.e. estimate using the whole spectrum)
262 casa::Float itsNoiseBox;
263
264 // true if median statistics is used as opposed to average of
265 // the lowest 80% of deviations (default)
266 casa::Bool itsUseMedian;
267
268 // true if spectra and mask data should be provided from
269 // scantable (default = true)
270 bool useScantable;
271
272 // shared object for nominal throw
273 casa::AipsError err ;
274};
275
276//
277///////////////////////////////////////////////////////////////////////////////
278
279} // namespace asap
280#endif // #ifndef STLINEFINDER_H
Note: See TracBrowser for help on using the repository browser.