casacore
FitGaussian.h
Go to the documentation of this file.
1//# FitGaussian.h: Multidimensional fitter class for Gaussians
2//# Copyright (C) 2001,2002
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$
28#ifndef SCIMATH_FITGAUSSIAN_H
29#define SCIMATH_FITGAUSSIAN_H
30
31#include <casacore/casa/aips.h>
32#include <casacore/casa/Arrays/Matrix.h>
33#include <casacore/casa/Logging/LogIO.h>
34
35namespace casacore { //# NAMESPACE CASACORE - BEGIN
36
37// <summary>Multidimensional fitter class for Gaussians.</summary>
38
39// <reviewed reviewer="" date="" tests="tFitGaussian">
40// </reviewed>
41
42// <prerequisite>
43// <li> <linkto class="Gaussian1D">Gaussian1D</linkto> class
44// <li> <linkto class="Gaussian2D">Gaussian2D</linkto> class
45// <li> <linkto class="Gaussian3D">Gaussian3D</linkto> class
46// <li> <linkto class="NonLinearFitLM">NonLinearFitLM</linkto> class
47// </prerequisite>
48
49// <etymology>
50// Fits Gaussians to data.
51// </etymology>
52
53// <synopsis>
54
55// <src>FitGaussian</src> is specially designed for fitting procedures in
56// code that must be generalized for general dimensionality and
57// number of components, and for complicated fits where the failure rate of
58// the standard nonlinear fitter is unacceptibly high.
59
60// <src>FitGaussian</src> essentially provides a Gaussian-adapted
61// interface for NonLinearFitLM. The user specifies the dimension,
62// number of gaussians, initial estimate, retry factors, and the data,
63// and the fitting proceeds automatically. Upon failure of the fitter it will
64// retry the fit according to the retry factors until a fit is completed
65// successfully. The user can optionally require as a criterion for success
66// that the RMS of the fit residuals not exceed some maximum value.
67
68// The retry factors are applied in different ways: the height and widths
69// are multiplied by the retry factors while the center and angles are
70// increased by their factors. As of 2002/07/12 these are applied randomly
71// (instead of sequentially) to different components and combinations of
72// components. The factors can be specified by the user, but a default
73// set is available. This random method is better than the sequential method
74// for a limited number of retries, but true optimization of the retry system
75// would demand the use of a more sophisticated method.
76// </synopsis>
77
78
79// <example>
80// <srcblock>
81// FitGaussian<Double> fitgauss(1,1);
82// Matrix<Double> x(5,1); x(0,0) = 0; x(1,0) = 1; x(2,0) = 2; x(3,0) = 3; x(4,0) = 4;
83// Vector<Double> y(5); y(0) = 0; y(1) = 1; y(2) = 4; y(3) = 1; y(4) = 1;
84// Matrix<Double> estimate(1,3);
85// estimate(0,0) = 1; estimate(0,1) = 1; estimate(0,2) = 1;
86// fitgauss.setFirstEstimate(estimate);
87// Matrix<Double> solution;
88// solution = fitgauss.fit(x,y);
89// cout << solution;
90// </srcblock>
91// </example>
92
93// <motivation>
94// Fitting multiple Gaussians is required for many different applications,
95// but requires a substantial amount of coding - especially if the
96// dimensionality of the image is not known to the programmer. Furthermore,
97// fitting multiple Gaussians has a very high failure rate. So, a specialized
98// Gaussian fitting class that retries from different initial estimates
99// until an acceptible fit was found was needed.
100// </motivation>
101
102// <templating arg=T>
103// <li> T must be a real data type compatible with NonLinearFitLM - Float or
104// Double.
105// </templating>
106
107// <thrown>
108// <li> AipsError if dimension is not 1, 2, or 3
109// <li> AipsError if incorrect parameter number specified.
110// <li> AipsError if estimate/retry/data arrays are of wrong dimension
111// </thrown>
112
113// <todo asof="2002/07/22">
114// <li> Optimize the default retry matrix
115// <li> Send fitting messages to logger instead of console
116// <li> Consider using a more sophisticated retry ststem (above).
117// <li> Check the estimates for reasonability, especially on failure of fit.
118// <li> Consider adding other models (polynomial, etc) to make this a Fit3D
119// class.
120// </todo>
121
122
123
124template <class T>
126{
127 public:
128
129 // Create the fitter. The dimension and the number of gaussians to fit
130 // can be modified later if necessary.
131 // <group>
133 FitGaussian(uInt dimension);
134 FitGaussian(uInt dimension, uInt numgaussians);
135 // </group>
136
137 // Adjust the number of dimensions
138 void setDimensions(uInt dimensions);
139
140 // Adjust the number of gaussians to fit
141 void setNumGaussians(uInt numgaussians);
142
143 // Set the initial estimate (the starting point of the first fit.)
144 void setFirstEstimate(const Matrix<T>& estimate);
145
146 // Set the maximum number of retries.
147 void setMaxRetries(uInt nretries) {itsMaxRetries = nretries;};
148
149 // Set the maximum amount of time to spend (in seconds). If time runs out
150 // during a fit the process will still complete that fit.
151 void setMaxTime(Double maxtime) {itsMaxTime = maxtime;};
152
153 // Set the retry factors, the values that are added/multiplied with the
154 // first estimate on subsequent attempts if the first attempt fails.
155 // Using the function with no argument sets the retry factors to the default.
156 // <group>
158 void setRetryFactors(const Matrix<T>& retryfactors);
159 // </group>
160
161 // Return the number of retry options available
163
164 // Mask out some parameters so that they are not modified during fitting
165 Bool &mask(uInt gaussian, uInt parameter);
166 const Bool &mask(uInt gaussian, uInt parameter) const;
167
168 // Run the fit, using the data provided in the arguments pos and f.
169 // The fit will retry from different initial estimates until it converges
170 // to a value with an RMS error less than maximumRMS. If this cannot be
171 // accomplished it will simply take the result that generated the best RMS.
172 Matrix<T> fit(const Matrix<T>& pos, const Vector<T>& f,
173 T maximumRMS = 1.0, uInt maxiter = 1024,
174 T convcriteria = 0.0001);
175 Matrix<T> fit(const Matrix<T>& pos,const Vector<T>& f,
176 const Vector<T>& sigma,
177 T maximumRMS = 1.0, uInt maxiter = 1024,
178 T convcriteria = 0.0001);
179
180 // Allow access to the fit parameters from this class
183
184 // Internal function for ensuring that parameters stay within their stated
185 // domains (see <src>Gaussian2D</src> and <src>Gaussian3D</src>.)
186 void correctParameters(Matrix<T>& parameters);
187
188 // Return the chi squared of the fit
190
191 // Return the RMS of the fit
192 T RMS();
193
194 // Returns True if the fit (eventually) converged to a value.
196
197
198 private:
199 uInt itsDimension; // how many dimensions (1, 2, or 3)
200 uInt itsNGaussians; // number of gaussians to fit
201 uInt itsMaxRetries; // maximum number of retries to attempt
202 Double itsMaxTime; // maximum time to spend fitting in secs
203 T itsChisquare; // chisquare of fit
204 T itsRMS; // RMS of fit (sqrt[chisquare / N])
205 Bool itsSuccess; // flags success or failure
207
208 Matrix<T> itsFirstEstimate; // user's estimate.
209 Matrix<T> itsRetryFctr; // source of retry information
210 Matrix<Bool> itsMask; // masks parameters not to change in fitting
211
212
213 // Sets the retry matrix to a default value. This is done automatically if
214 // the retry matrix is not set directly.
216
217 //Add one or more rows to the retry matrix.
218 void expandRetryMatrix(uInt rowstoadd);
219
220 //Find the number of unmasked parameters to be fit
222
223 // The solutions to the fit
225
226 // The errors on the solution parameters
228
229};
230
231
232
233} //# NAMESPACE CASACORE - END
234
235#ifndef CASACORE_NO_AUTO_TEMPLATES
236#include <casacore/scimath/Fitting/FitGaussian.tcc>
237#endif //# CASACORE_NO_AUTO_TEMPLATES
238#endif
239
240
241
242
243
244
245
246
Matrix< T > fit(const Matrix< T > &pos, const Vector< T > &f, T maximumRMS=1.0, uInt maxiter=1024, T convcriteria=0.0001)
Run the fit, using the data provided in the arguments pos and f.
Bool converged()
Returns True if the fit (eventually) converged to a value.
uInt nRetryFactors()
Return the number of retry options available.
Definition: FitGaussian.h:162
const Matrix< T > & solution()
Allow access to the fit parameters from this class.
Definition: FitGaussian.h:181
T chisquared()
Return the chi squared of the fit.
Bool & mask(uInt gaussian, uInt parameter)
Mask out some parameters so that they are not modified during fitting.
Matrix< T > itsSolutionErrors
The errors on the solution parameters.
Definition: FitGaussian.h:227
const Matrix< T > & errors()
Definition: FitGaussian.h:182
void setRetryFactors(const Matrix< T > &retryfactors)
Matrix< Bool > itsMask
Definition: FitGaussian.h:210
void expandRetryMatrix(uInt rowstoadd)
Add one or more rows to the retry matrix.
uInt countFreeParameters()
Find the number of unmasked parameters to be fit.
FitGaussian(uInt dimension, uInt numgaussians)
T RMS()
Return the RMS of the fit.
Matrix< T > itsFirstEstimate
Definition: FitGaussian.h:208
Matrix< T > defaultRetryMatrix()
Sets the retry matrix to a default value.
Matrix< T > itsRetryFctr
Definition: FitGaussian.h:209
void setFirstEstimate(const Matrix< T > &estimate)
Set the initial estimate (the starting point of the first fit.)
FitGaussian(uInt dimension)
void setMaxTime(Double maxtime)
Set the maximum amount of time to spend (in seconds).
Definition: FitGaussian.h:151
void setMaxRetries(uInt nretries)
Set the maximum number of retries.
Definition: FitGaussian.h:147
const Bool & mask(uInt gaussian, uInt parameter) const
void correctParameters(Matrix< T > &parameters)
Internal function for ensuring that parameters stay within their stated domains (see Gaussian2D and G...
void setNumGaussians(uInt numgaussians)
Adjust the number of gaussians to fit.
FitGaussian()
Create the fitter.
Matrix< T > itsSolutionParameters
The solutions to the fit.
Definition: FitGaussian.h:224
void setDimensions(uInt dimensions)
Adjust the number of dimensions.
void setRetryFactors()
Set the retry factors, the values that are added/multiplied with the first estimate on subsequent att...
Matrix< T > fit(const Matrix< T > &pos, const Vector< T > &f, const Vector< T > &sigma, T maximumRMS=1.0, uInt maxiter=1024, T convcriteria=0.0001)
size_t nrow() const
The number of rows in the Matrix, i.e.
Definition: Matrix.h:300
this file contains all the compiler specific defines
Definition: mainpage.dox:28
unsigned int uInt
Definition: aipstype.h:51
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
double Double
Definition: aipstype.h:55