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