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 
35 namespace 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 
124 template <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
182  const Matrix<T> &errors(){return itsSolutionErrors;};
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 
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
T chisquared()
Return the chi squared of the fit.
const Matrix< T > & solution()
Allow access to the fit parameters from this class.
Definition: FitGaussian.h:181
Matrix< T > itsSolutionErrors
The errors on the solution parameters.
Definition: FitGaussian.h:227
void setRetryFactors(const Matrix< T > &retryfactors)
Matrix< T > defaultRetryMatrix()
Sets the retry matrix to a default value.
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.
Bool & mask(uInt gaussian, uInt parameter)
Mask out some parameters so that they are not modified during fitting.
FitGaussian(uInt dimension, uInt numgaussians)
T RMS()
Return the RMS of the fit.
Matrix< T > itsFirstEstimate
Definition: FitGaussian.h:208
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
const Bool & mask(uInt gaussian, uInt parameter) const
void setMaxRetries(uInt nretries)
Set the maximum number of retries.
Definition: FitGaussian.h:147
void correctParameters(Matrix< T > &parameters)
Internal function for ensuring that parameters stay within their stated domains (see Gaussian2D and G...
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)
const Matrix< T > & errors()
Definition: FitGaussian.h:182
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.
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...
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