source: trunk/src/STLineFinder.h@ 2046

Last change on this file since 2046 was 2012, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



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