source: branches/casa-release-4_3-test/external-alma/components/SpectralComponents/SpectralEstimate.h@ 3040

Last change on this file since 3040 was 2980, checked in by Malte Marquarding, 10 years ago

Add a copy of casacore/components/SpectralComponents to external-alma directory to prepare for its removal from casacore-trunk. DOn;t activate in SConscript yet.

File size: 7.8 KB
Line 
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
36namespace casa { //# NAMESPACE CASA - BEGIN
37
38//# Forward Declarations
39class GaussianSpectralElement;
40template <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
95class 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
Note: See TracBrowser for help on using the repository browser.