FitGaussian.h

Go to the documentation of this file.

// # FitGaussian.h: Multidimensional fitter class for Gaussians
// # Copyright (C) 2001,2002
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library is distributed in the hope that it will be useful, but WITHOUT
// # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
#ifndef SCIMATH_FITGAUSSIAN_H
#define SCIMATH_FITGAUSSIAN_H
 
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/casa/Logging/LogIO.h>
 
namespace casacore {  // # NAMESPACE CASACORE - BEGIN
 
// <summary>Multidimensional fitter class for Gaussians.</summary>
 
// <reviewed reviewer="" date="" tests="tFitGaussian">
// </reviewed>
 
// <prerequisite>
//   <li> <linkto class="Gaussian1D">Gaussian1D</linkto> class
//   <li> <linkto class="Gaussian2D">Gaussian2D</linkto> class
//   <li> <linkto class="Gaussian3D">Gaussian3D</linkto> class
//   <li> <linkto class="NonLinearFitLM">NonLinearFitLM</linkto> class
// </prerequisite>
 
// <etymology>
// Fits Gaussians to data.
// </etymology>
 
// <synopsis>
 
// <src>FitGaussian</src> is specially designed for fitting procedures in
// code that must be generalized for general dimensionality and
// number of components, and for complicated fits where the failure rate of
// the standard nonlinear fitter is unacceptibly high.
 
// <src>FitGaussian</src> essentially provides a Gaussian-adapted
// interface for NonLinearFitLM.  The user specifies the dimension,
// number of gaussians, initial estimate, retry factors, and the data,
// and the fitting proceeds automatically.  Upon failure of the fitter it will
// retry the fit according to the retry factors until a fit is completed
// successfully.  The user can optionally require as a criterion for success
// that the RMS of the fit residuals not exceed some maximum value.
 
// The retry factors are applied in different ways: the height and widths
// are multiplied by the retry factors while the center and angles are
// increased by their factors.  As of 2002/07/12 these are applied randomly
// (instead of sequentially) to different components and combinations of
// components.  The factors can be specified by the user, but a default
// set is available.  This random method is better than the sequential method
// for a limited number of retries, but true optimization of the retry system
// would demand the use of a more sophisticated method.
// </synopsis>
 
// <example>
// <srcblock>
// FitGaussian<Double> fitgauss(1,1);
// 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;
// Vector<Double> y(5); y(0) = 0; y(1) = 1; y(2) = 4; y(3) = 1; y(4) = 1;
// Matrix<Double> estimate(1,3);
// estimate(0,0) = 1; estimate(0,1) = 1; estimate(0,2) = 1;
// fitgauss.setFirstEstimate(estimate);
// Matrix<Double> solution;
// solution = fitgauss.fit(x,y);
// cout << solution;
// </srcblock>
// </example>
 
// <motivation>
// Fitting multiple Gaussians is required for many different applications,
// but requires a substantial amount of coding - especially if the
// dimensionality of the image is not known to the programmer.  Furthermore,
// fitting multiple Gaussians has a very high failure rate.  So, a specialized
// Gaussian fitting class that retries from different initial estimates
// until an acceptible fit was found was needed.
// </motivation>
 
// <templating arg=T>
//  <li> T must be a real data type compatible with NonLinearFitLM - Float or
//  Double.
// </templating>
 
// <thrown>
//    <li> AipsError if dimension is not 1, 2, or 3
//    <li> AipsError if incorrect parameter number specified.
//    <li> AipsError if estimate/retry/data arrays are of wrong dimension
// </thrown>
 
// <todo asof="2002/07/22">
//   <li> Optimize the default retry matrix
//   <li> Send fitting messages to logger instead of console
//   <li> Consider using a more sophisticated retry ststem (above).
//   <li> Check the estimates for reasonability, especially on failure of fit.
//   <li> Consider adding other models (polynomial, etc) to make this a Fit3D
//        class.
// </todo>
 
template <class T>
class FitGaussian {
 public:
  // Create the fitter.  The dimension and the number of gaussians to fit
  // can be modified later if necessary.
  // <group>
  FitGaussian();
  FitGaussian(uInt dimension);
  FitGaussian(uInt dimension, uInt numgaussians);
  // </group>
 
  // Adjust the number of dimensions
  void setDimensions(uInt dimensions);
 
  // Adjust the number of gaussians to fit
  void setNumGaussians(uInt numgaussians);
 
  // Set the initial estimate (the starting point of the first fit.)
  void setFirstEstimate(const Matrix<T>& estimate);
 
  // Set the maximum number of retries.
  void setMaxRetries(uInt nretries) { itsMaxRetries = nretries; };
 
  // Set the maximum amount of time to spend (in seconds).  If time runs out
  // during a fit the process will still complete that fit.
  void setMaxTime(Double maxtime) { itsMaxTime = maxtime; };
 
  // Set the retry factors, the values that are added/multiplied with the
  // first estimate on subsequent attempts if the first attempt fails.
  // Using the function with no argument sets the retry factors to the default.
  // <group>
  void setRetryFactors();
  void setRetryFactors(const Matrix<T>& retryfactors);
  // </group>
 
  // Return the number of retry options available
  uInt nRetryFactors() { return itsRetryFctr.nrow(); };
 
  // Mask out some parameters so that they are not modified during fitting
  Bool& mask(uInt gaussian, uInt parameter);
  const Bool& mask(uInt gaussian, uInt parameter) const;
 
  // Run the fit, using the data provided in the arguments pos and f.
  // The fit will retry from different initial estimates until it converges
  // to a value with an RMS error less than maximumRMS.  If this cannot be
  // accomplished it will simply take the result that generated the best RMS.
  Matrix<T> fit(const Matrix<T>& pos, const Vector<T>& f, T maximumRMS = 1.0, uInt maxiter = 1024,
                T convcriteria = 0.0001);
  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);
 
  // Allow access to the fit parameters from this class
  const Matrix<T>& solution() { return itsSolutionParameters; };
  const Matrix<T>& errors() { return itsSolutionErrors; };
 
  // Internal function for ensuring that parameters stay within their stated
  // domains (see <src>Gaussian2D</src> and <src>Gaussian3D</src>.)
  void correctParameters(Matrix<T>& parameters);
 
  // Return the chi squared of the fit
  T chisquared();
 
  // Return the RMS of the fit
  T RMS();
 
  // Returns True if the fit (eventually) converged to a value.
  Bool converged();
 
 private:
  uInt itsDimension;   // how many dimensions (1, 2, or 3)
  uInt itsNGaussians;  // number of gaussians to fit
  uInt itsMaxRetries;  // maximum number of retries to attempt
  Double itsMaxTime;   // maximum time to spend fitting in secs
  T itsChisquare;      // chisquare of fit
  T itsRMS;            // RMS of fit (sqrt[chisquare / N])
  Bool itsSuccess;     // flags success or failure
  LogIO os;
 
  Matrix<T> itsFirstEstimate;  // user's estimate.
  Matrix<T> itsRetryFctr;      // source of retry information
  Matrix<Bool> itsMask;        // masks parameters not to change in fitting
 
  // Sets the retry matrix to a default value.  This is done automatically if
  // the retry matrix is not set directly.
  Matrix<T> defaultRetryMatrix();
 
  // Add one or more rows to the retry matrix.
  void expandRetryMatrix(uInt rowstoadd);
 
  // Find the number of unmasked parameters to be fit
  uInt countFreeParameters();
 
  // The solutions to the fit
  Matrix<T> itsSolutionParameters;
 
  // The errors on the solution parameters
  Matrix<T> itsSolutionErrors;
};
 
}  // namespace casacore
 
#ifndef CASACORE_NO_AUTO_TEMPLATES
#include <casacore/scimath/Fitting/FitGaussian.tcc>
#endif  // # CASACORE_NO_AUTO_TEMPLATES
#endif