casacore
Loading...
Searching...
No Matches
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.
const Matrix< T > & solution()
Allow access to the fit parameters from this class.
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.
const Matrix< T > & errors()
void setRetryFactors(const Matrix< T > &retryfactors)
Matrix< Bool > itsMask
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
Matrix< T > defaultRetryMatrix()
Sets the retry matrix to a default value.
Matrix< T > itsRetryFctr
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).
void setMaxRetries(uInt nretries)
Set the maximum number of retries.
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.
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