1 | //# SpectralEstimate.h: Get an initial estimate for spectral lines
|
---|
2 | //# Copyright (C) 2001,2002,2003,2004
|
---|
3 | //# Associated Universities, Inc. Washington DC, USA.
|
---|
4 | //#
|
---|
5 | //# This library is free software; you can redistribute it and/or modify it
|
---|
6 | //# under the terms of the GNU Library General Public License as published by
|
---|
7 | //# the Free Software Foundation; either version 2 of the License, or (at your
|
---|
8 | //# option) any later version.
|
---|
9 | //#
|
---|
10 | //# This library is distributed in the hope that it will be useful, but WITHOUT
|
---|
11 | //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
---|
12 | //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
|
---|
13 | //# License for more details.
|
---|
14 | //#
|
---|
15 | //# You should have received a copy of the GNU Library General Public License
|
---|
16 | //# along with this library; if not, write to the Free Software Foundation,
|
---|
17 | //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
|
---|
18 | //#
|
---|
19 | //# Correspondence concerning AIPS++ should be addressed as follows:
|
---|
20 | //# Internet email: aips2-request@nrao.edu.
|
---|
21 | //# Postal address: AIPS++ Project Office
|
---|
22 | //# National Radio Astronomy Observatory
|
---|
23 | //# 520 Edgemont Road
|
---|
24 | //# Charlottesville, VA 22903-2475 USA
|
---|
25 | //#
|
---|
26 | //#
|
---|
27 | //# $Id: SpectralEstimate.h 21229 2012-04-02 12:00:20Z gervandiepen $
|
---|
28 |
|
---|
29 | #ifndef COMPONENTS_SPECTRALESTIMATE_H
|
---|
30 | #define COMPONENTS_SPECTRALESTIMATE_H
|
---|
31 |
|
---|
32 | #include <casa/aips.h>
|
---|
33 | #include <components/SpectralComponents/SpectralElement.h>
|
---|
34 | #include <components/SpectralComponents/SpectralList.h>
|
---|
35 |
|
---|
36 | namespace casa { //# NAMESPACE CASA - BEGIN
|
---|
37 |
|
---|
38 | //# Forward Declarations
|
---|
39 | class GaussianSpectralElement;
|
---|
40 | template <class T> class Vector;
|
---|
41 |
|
---|
42 | // <summary>
|
---|
43 | // Get an initial estimate for spectral lines
|
---|
44 | // </summary>
|
---|
45 |
|
---|
46 | // <use visibility=export>
|
---|
47 |
|
---|
48 | // <reviewed reviewer="" date="yyyy/mm/dd" tests="tSpectralFit" demos="">
|
---|
49 | // </reviewed>
|
---|
50 |
|
---|
51 | // <prerequisite>
|
---|
52 | // <li> <linkto class=SpectralFit>SpectralFit</linkto> class
|
---|
53 | // </prerequisite>
|
---|
54 | //
|
---|
55 | // <etymology>
|
---|
56 | // From spectral line and estimate
|
---|
57 | // </etymology>
|
---|
58 | //
|
---|
59 | // <synopsis>
|
---|
60 | // The SpectralEstimate class obtains an initial guess for spectral
|
---|
61 | // components. The current implementation uses the entire
|
---|
62 | // profile as signal region, or can set a window to be searched around
|
---|
63 | // the highest peak automatically. A window can also be set manually.
|
---|
64 | // The second derivative of
|
---|
65 | // the profile in the signal region is calculated by fitting
|
---|
66 | // a second degree polynomal. The smoothing parameter Q
|
---|
67 | // determines the number of points used for this (=2*Q+1).
|
---|
68 | // The gaussians can then be estimated as described by
|
---|
69 | // Schwarz, 1968, Bull.Astr.Inst.Netherlands, Volume 19, 405.
|
---|
70 | //
|
---|
71 | // The elements guessed can be used in the
|
---|
72 | // <linkto class=SpectralFit>SpectralFit</linkto> class.
|
---|
73 | //
|
---|
74 | // The default type found is a Gaussian, defined as:
|
---|
75 | // <srcblock>
|
---|
76 | // AMPL.exp[ -(x-CENTER)<sup>2</sup>/2 SIGMA<sup>2</sup>]
|
---|
77 | // </srcblock>
|
---|
78 | //
|
---|
79 | // The parameter estimates are returned in units of zero-based
|
---|
80 | // pixel indices.
|
---|
81 | // </synopsis>
|
---|
82 | //
|
---|
83 | // <example>
|
---|
84 | // </example>
|
---|
85 | //
|
---|
86 | // <motivation>
|
---|
87 | // To have an automatic method to find spectral lines
|
---|
88 | // </motivation>
|
---|
89 | //
|
---|
90 | // <todo asof="2001/02/14">
|
---|
91 | // <li> find a way to get to absorption lines as well
|
---|
92 | // <li> add more estimation options
|
---|
93 | // </todo>
|
---|
94 |
|
---|
95 | class SpectralEstimate {
|
---|
96 | public:
|
---|
97 | //# Constants
|
---|
98 | // Default maximum number of components to be found
|
---|
99 | static const uInt MAXPAR = 200;
|
---|
100 | //# Enumerations
|
---|
101 | //# Friends
|
---|
102 |
|
---|
103 | //# Constructors
|
---|
104 | // Default constructor creates a default estimator (default max number
|
---|
105 | // of components to be found is 200) with the given maximum number
|
---|
106 | // of components that will be found. A value of zero will indicate
|
---|
107 | // an unlimited number.
|
---|
108 | explicit SpectralEstimate(const uInt maxpar=MAXPAR);
|
---|
109 | // Create an estimator with the given maximum number of possible
|
---|
110 | // elements. A value of zero will indicate an unlimited number.
|
---|
111 | // Construct with a given rms in profiles, a cutoff for amplitudes
|
---|
112 | // found, and a minimum width. Cutoff and minsigma default to 0.0, maximum
|
---|
113 | // size of list produced to 200.
|
---|
114 | explicit SpectralEstimate(const Double rms,
|
---|
115 | const Double cutoff=0.0, const Double minsigma=0.0,
|
---|
116 | const uInt maxpar=MAXPAR);
|
---|
117 | // Copy constructor (deep copy)
|
---|
118 | SpectralEstimate(const SpectralEstimate &other);
|
---|
119 |
|
---|
120 | //#Destructor
|
---|
121 | // Destructor
|
---|
122 | ~SpectralEstimate();
|
---|
123 |
|
---|
124 | //# Operators
|
---|
125 | // Assignment (copy semantics)
|
---|
126 | SpectralEstimate &operator=(const SpectralEstimate &other);
|
---|
127 |
|
---|
128 | //# Member functions
|
---|
129 | // Generate the estimates for a profile and return the
|
---|
130 | // list found. The first function returns component parameters
|
---|
131 | // in units of pixel indices. The second function calls the first
|
---|
132 | // and then converts to the specified abcissa space (the supplied
|
---|
133 | // vector must be monotonic); if the pixel-based center is out of range
|
---|
134 | // of the supplied abcissa vector the conversion is done via extrapolation.
|
---|
135 | // The der pointer is meant for debugging, and can return
|
---|
136 | // the derivative profile. The second function throws an AipsError
|
---|
137 | // if the vectors are not the same length.
|
---|
138 | // <group>
|
---|
139 | template <class MT>
|
---|
140 | const SpectralList& estimate(const Vector<MT>& ordinate,
|
---|
141 | Vector<MT> *der = 0);
|
---|
142 | template <class MT>
|
---|
143 | const SpectralList& estimate(const Vector<MT>& abcissa,
|
---|
144 | const Vector<MT>& ordinate);
|
---|
145 | // </group>
|
---|
146 |
|
---|
147 | // Return the list found.
|
---|
148 | const SpectralList &list() const {return slist_p; };
|
---|
149 |
|
---|
150 | // Set estimation parameters
|
---|
151 | // <group>
|
---|
152 | // Set the profile's estimated rms (forced to abs(rms))
|
---|
153 | void setRMS(const Double rms=0.0);
|
---|
154 | // Set the amplitude cutoff for valid estimate (forced to max(0,cutoff))
|
---|
155 | void setCutoff(const Double cutoff=0.0);
|
---|
156 | // Set the minimum width allowed (forced to max(0,minsigma))
|
---|
157 | void setMinSigma(const Double minsigma=0.0);
|
---|
158 | // Set the number of points consider at each side of test point (i.e. a
|
---|
159 | // width of 2q+1 is taken). Default internally is 2; max(1,q) taken.
|
---|
160 | void setQ(const uInt q=2);
|
---|
161 | // Set a region [lo,hi] over which to estimate. Lo and hi are given as
|
---|
162 | // zero-based vector indices.
|
---|
163 | void setRegion(const Int lo, const Int hi);
|
---|
164 | // Do you want to look in an automatically determined window with signal?
|
---|
165 | // Default is False, meaning the full (possibly regioned) profile.
|
---|
166 | void setWindowing(const Bool win=False);
|
---|
167 | // Set the maximum number of estimates to find (forced to >=1; 200 default)
|
---|
168 | void setMaxN(const uInt maxpar=MAXPAR);
|
---|
169 | // </group>
|
---|
170 |
|
---|
171 | private:
|
---|
172 | //#Data
|
---|
173 | // Use window search
|
---|
174 | Bool useWindow_p;
|
---|
175 | // rms estimate in profile
|
---|
176 | Double rms_p;
|
---|
177 | // Source cutoff amplitude
|
---|
178 | Double cutoff_p;
|
---|
179 | // Window low and end value
|
---|
180 | // <group>
|
---|
181 | Int windowLow_p;
|
---|
182 | Int windowEnd_p;
|
---|
183 | // </group>
|
---|
184 | // Region low and high value
|
---|
185 | // <group>
|
---|
186 | Int regionLow_p;
|
---|
187 | Int regionEnd_p;
|
---|
188 | // </group>
|
---|
189 | // Smoothing parameter. I.e. 2q+1 points are taken
|
---|
190 | Int q_p;
|
---|
191 | // Internal cashing of calculated values based on q
|
---|
192 | // <group>
|
---|
193 | Double a_p;
|
---|
194 | Double b_p;
|
---|
195 | // </group>
|
---|
196 | // The minimum Gaussian width
|
---|
197 | Double sigmin_p;
|
---|
198 | // The second derivatives
|
---|
199 | Double *deriv_p;
|
---|
200 | // The list of components
|
---|
201 | SpectralList slist_p;
|
---|
202 | // The length of the current profile being estimated
|
---|
203 | uInt lprof_p;
|
---|
204 |
|
---|
205 | //# Member functions
|
---|
206 | // Get the window or the total spectrum
|
---|
207 | template <class MT>
|
---|
208 | uInt window(const Vector<MT> &prof);
|
---|
209 | // Get the second derivatives
|
---|
210 | template <class MT>
|
---|
211 | void findc2(const Vector<MT> &prof);
|
---|
212 | // Find the Gaussians
|
---|
213 | template <class MT>
|
---|
214 | void findga(const Vector<MT> &prof);
|
---|
215 | // Convert the parameters of the components in the list from
|
---|
216 | // pixel-based indices to the given abcissa-vector space.
|
---|
217 | template <class MT> GaussianSpectralElement convertElement (const Vector<MT>& abcissa,
|
---|
218 | const GaussianSpectralElement& el) const;
|
---|
219 | };
|
---|
220 |
|
---|
221 |
|
---|
222 | } //# NAMESPACE CASA - END
|
---|
223 |
|
---|
224 | #ifndef CASACORE_NO_AUTO_TEMPLATES
|
---|
225 | #include <components/SpectralComponents/Spectral2Estimate.tcc>
|
---|
226 | #endif //# CASACORE_NO_AUTO_TEMPLATES
|
---|
227 | #endif
|
---|