casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
FitGaussian.h
Go to the documentation of this file.
00001 //# FitGaussian.h: Multidimensional fitter class for Gaussians
00002 //# Copyright (C) 2001,2002
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //#
00027 //# $Id: FitGaussian.h 20229 2008-01-29 15:19:06Z gervandiepen $
00028 #ifndef SCIMATH_FITGAUSSIAN_H
00029 #define SCIMATH_FITGAUSSIAN_H
00030 
00031 #include <casa/aips.h>
00032 #include <casa/Arrays/Matrix.h>
00033 #include <casa/Logging/LogIO.h>
00034 
00035 namespace casa { //# NAMESPACE CASA - BEGIN
00036 
00037 // <summary>Multidimensional fitter class for Gaussians.</summary>
00038 
00039 // <reviewed reviewer="" date="" tests="tFitGaussian">
00040 // </reviewed>
00041 
00042 // <prerequisite>
00043 //   <li> <linkto class="Gaussian1D">Gaussian1D</linkto> class
00044 //   <li> <linkto class="Gaussian2D">Gaussian2D</linkto> class
00045 //   <li> <linkto class="Gaussian3D">Gaussian3D</linkto> class
00046 //   <li> <linkto class="NonLinearFitLM">NonLinearFitLM</linkto> class
00047 // </prerequisite>
00048 
00049 // <etymology>
00050 // Fits Gaussians to data.
00051 // </etymology>
00052 
00053 // <synopsis>
00054 
00055 // <src>FitGaussian</src> is specially designed for fitting procedures in
00056 // code that must be generalized for general dimensionality and 
00057 // number of components, and for complicated fits where the failure rate of
00058 // the standard nonlinear fitter is unacceptibly high.
00059 
00060 // <src>FitGaussian</src> essentially provides a Gaussian-adapted 
00061 // interface for NonLinearFitLM.  The user specifies the dimension, 
00062 // number of gaussians, initial estimate, retry factors, and the data,
00063 // and the fitting proceeds automatically.  Upon failure of the fitter it will 
00064 // retry the fit according to the retry factors until a fit is completed
00065 // successfully.  The user can optionally require as a criterion for success
00066 // that the RMS of the fit residuals not exceed some maximum value.
00067 
00068 // The retry factors are applied in different ways: the height and widths
00069 // are multiplied by the retry factors while the center and angles are
00070 // increased by their factors.  As of 2002/07/12 these are applied randomly
00071 // (instead of sequentially) to different components and combinations of
00072 // components.  The factors can be specified by the user, but a default
00073 // set is available.  This random method is better than the sequential method
00074 // for a limited number of retries, but true optimization of the retry system
00075 // would demand the use of a more sophisticated method.
00076 // </synopsis>
00077 
00078 
00079 // <example>
00080 // <srcblock>
00081 // FitGaussian<Double> fitgauss(1,1);
00082 // 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;
00083 // Vector<Double> y(5); y(0) = 0; y(1) = 1; y(2) = 4; y(3) = 1; y(4) = 1;
00084 // Matrix<Double> estimate(1,3);
00085 // estimate(0,0) = 1; estimate(0,1) = 1; estimate(0,2) = 1;
00086 // fitgauss.setFirstEstimate(estimate);
00087 // Matrix<Double> solution;
00088 // solution = fitgauss.fit(x,y);
00089 // cout << solution;
00090 // </srcblock>
00091 // </example>
00092 
00093 // <motivation>
00094 // Fitting multiple Gaussians is required for many different applications,
00095 // but requires a substantial amount of coding - especially if the 
00096 // dimensionality of the image is not known to the programmer.  Furthermore,
00097 // fitting multiple Gaussians has a very high failure rate.  So, a specialized
00098 // Gaussian fitting class that retries from different initial estimates
00099 // until an acceptible fit was found was needed.
00100 // </motivation>
00101 
00102 // <templating arg=T>
00103 //  <li> T must be a real data type compatible with NonLinearFitLM - Float or
00104 //  Double.
00105 // </templating>
00106 
00107 // <thrown>
00108 //    <li> AipsError if dimension is not 1, 2, or 3
00109 //    <li> AipsError if incorrect parameter number specified.
00110 //    <li> AipsError if estimate/retry/data arrays are of wrong dimension
00111 // </thrown>
00112   
00113 // <todo asof="2002/07/22">
00114 //   <li> Optimize the default retry matrix
00115 //   <li> Send fitting messages to logger instead of console
00116 //   <li> Consider using a more sophisticated retry ststem (above).
00117 //   <li> Check the estimates for reasonability, especially on failure of fit.
00118 //   <li> Consider adding other models (polynomial, etc) to make this a Fit3D
00119 //        class.
00120 // </todo>
00121 
00122 
00123 
00124 template <class T>
00125 class FitGaussian
00126 {
00127   public:
00128 
00129   // Create the fitter.  The dimension and the number of gaussians to fit
00130   // can be modified later if necessary.
00131   // <group>
00132   FitGaussian();
00133   FitGaussian(uInt dimension);
00134   FitGaussian(uInt dimension, uInt numgaussians);
00135   // </group>
00136 
00137   // Adjust the number of dimensions
00138   void setDimensions(uInt dimensions);
00139 
00140   // Adjust the number of gaussians to fit
00141   void setNumGaussians(uInt numgaussians);
00142 
00143   // Set the initial estimate (the starting point of the first fit.)
00144   void setFirstEstimate(const Matrix<T>& estimate);
00145 
00146   // Set the maximum number of retries.
00147   void setMaxRetries(uInt nretries) {itsMaxRetries = nretries;};
00148 
00149   // Set the maximum amount of time to spend (in seconds).  If time runs out
00150   // during a fit the process will still complete that fit.
00151   void setMaxTime(Double maxtime) {itsMaxTime = maxtime;};
00152 
00153   // Set the retry factors, the values that are added/multiplied with the
00154   // first estimate on subsequent attempts if the first attempt fails.
00155   // Using the function with no argument sets the retry factors to the default.
00156   // <group>
00157   void setRetryFactors();
00158   void setRetryFactors(const Matrix<T>& retryfactors);
00159   // </group>
00160 
00161   // Return the number of retry options available
00162   uInt nRetryFactors() {return itsRetryFctr.nrow();};
00163 
00164   // Mask out some parameters so that they are not modified during fitting
00165   Bool &mask(uInt gaussian, uInt parameter);
00166   const Bool &mask(uInt gaussian, uInt parameter) const;
00167 
00168   // Run the fit, using the data provided in the arguments pos and f.
00169   // The fit will retry from different initial estimates until it converges
00170   // to a value with an RMS error less than maximumRMS.  If this cannot be
00171   // accomplished it will simply take the result that generated the best RMS.
00172   Matrix<T> fit(const Matrix<T>& pos, const Vector<T>& f,
00173                 T maximumRMS = 1.0, uInt maxiter = 1024, 
00174                 T convcriteria = 0.0001);
00175   Matrix<T> fit(const Matrix<T>& pos,const Vector<T>& f,
00176                 const Vector<T>& sigma,
00177                 T maximumRMS = 1.0, uInt maxiter = 1024, 
00178                 T convcriteria = 0.0001);
00179 
00180   // Internal function for ensuring that parameters stay within their stated
00181   // domains (see <src>Gaussian2D</src> and <src>Gaussian3D</src>.)
00182   void correctParameters(Matrix<T>& parameters);
00183 
00184   // Return the chi squared of the fit
00185   T chisquared();
00186 
00187   // Return the RMS of the fit
00188   T RMS();
00189 
00190   // Returns True if the fit (eventually) converged to a value.
00191   Bool converged();
00192 
00193 
00194   private:
00195   uInt itsDimension;           // how many dimensions (1, 2, or 3)
00196   uInt itsNGaussians;          // number of gaussians to fit
00197   uInt itsMaxRetries;          // maximum number of retries to attempt
00198   Double itsMaxTime;           // maximum time to spend fitting in secs
00199   T itsChisquare;              // chisquare of fit
00200   T itsRMS;                    // RMS of fit (sqrt[chisquare / N])
00201   Bool itsSuccess;             // flags success or failure
00202   LogIO os;
00203 
00204   Matrix<T> itsFirstEstimate;  // user's estimate.
00205   Matrix<T> itsRetryFctr;      // source of retry information
00206   Matrix<Bool> itsMask;        // masks parameters not to change in fitting
00207 
00208   
00209   // Sets the retry matrix to a default value.  This is done automatically if
00210   // the retry matrix is not set directly.
00211   Matrix<T> defaultRetryMatrix();
00212 
00213   //Add one or more rows to the retry matrix.
00214   void expandRetryMatrix(uInt rowstoadd);
00215 
00216   //Find the number of unmasked parameters to be fit
00217   uInt countFreeParameters();
00218 };
00219 
00220 
00221 
00222 } //# NAMESPACE CASA - END
00223 
00224 #ifndef CASACORE_NO_AUTO_TEMPLATES
00225 #include <scimath/Fitting/FitGaussian.tcc>
00226 #endif //# CASACORE_NO_AUTO_TEMPLATES
00227 #endif
00228 
00229 
00230 
00231 
00232 
00233 
00234 
00235