Fit2D.h

Go to the documentation of this file.
00001 //# Fit2D.h: Class to fit 2-D objects to Lattices or Arrays
00002 //# Copyright (C) 1997,1998,1999,2000,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 //#   $Id$
00027 
00028 #ifndef LATTICES_FIT2D_H
00029 #define LATTICES_FIT2D_H
00030 
00031 //# Includes
00032 #include <casacore/casa/aips.h>
00033 #include <casacore/scimath/Functionals/CompoundFunction.h>
00034 #include <casacore/casa/BasicSL/Constants.h>
00035 #include <casacore/scimath/Fitting/NonLinearFitLM.h>
00036 #include <casacore/casa/Logging/LogIO.h>
00037 
00038 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00039 
00040 template<class T> class Array;
00041 template<class T> class Matrix;
00042 template<class T> class Vector;
00043 template<class T> class Lattice;
00044 template<class T> class MaskedLattice;
00045 
00046 
00047 // <summary>
00048 // Fit 2-D objects to 2-D Lattices or Arrays
00049 // </summary>
00050 
00051 // <use visibility=export>
00052 
00053 // <reviewed reviewer="" date="" tests="">
00054 // </reviewed>
00055 
00056 // <prerequisite>
00057 //   <li> <linkto class=Lattice>Lattice</linkto>
00058 // </prerequisite>
00059 
00060 // <synopsis> 
00061 // This class allows you to fit different types of 2-D models
00062 // to either Lattices or Arrays.  These must be 2 dimensional;
00063 // for Lattices, the appropriate 2-D Lattice can be made with
00064 // the SubLattice class.
00065 //
00066 // You may fit more than one model simultaneously to the data.
00067 // Models are added with the addModel method.   With this method,
00068 // you also specify the initial guesses of the parameters of
00069 // the model.    Any parameters involving coordinates are
00070 // expected in zero-relative absolute pixel coordinates (e.g. the centre of
00071 // a model).  Additionally with the addModel method, 
00072 // you may specify which parameters are to be held fixed
00073 // during the fitting process.  This is done with the 
00074 // parameterMask Vector which is in the same order as the
00075 // parameter Vector.  A value of True indicates the parameter
00076 // will be fitted for.  Presently, when you say fix the minor axis,
00077 // you really end up fixing the axial ratio (internals).  I don't
00078 // have a solution for this presently.
00079 // 
00080 // For Gaussians, the parameter Vector (input or output) consists, in order, of
00081 // the peak, x location, y location, FWHM of major axis, FWHM of minor axis, 
00082 // and position angle of the major axis (in radians). The 
00083 // position angle is positive +x to +y 
00084 // in the pixel coordinate system ([0,0] in center of image) and 
00085 // in the range -2pi to 2pi.  When the solution is recovered, the
00086 // position angle will be in the range 0 to pi.
00087 //
00088 // </synopsis> 
00089 // <example>
00090 // <srcblock>
00091 // </srcblock>
00092 // </example>
00093 
00094 // <todo asof="1998/12/11">
00095 //  <li> template it 
00096 //  <li> Speed up some Array calculations indexed with IPositions
00097 //  <li> Don't handle Lattices simply by getting pixels into Arrays
00098 //  <li> make an addModel interface taking functionals
00099 // </todo>
00100 
00101 class Fit2D 
00102 {
00103 public:
00104 
00105     // Enum describing the different models you can fit
00106     enum Types {
00107       GAUSSIAN = 0,
00108       DISK = 1,
00109       LEVEL=2,
00110       PLANE=3,
00111       nTypes
00112     };
00113 
00114     // Enum describing output error conditions
00115     enum ErrorTypes {
00116 // ok
00117       OK = 0,
00118 // Did not converge
00119       NOCONVERGE = 1,
00120 // Solution failed
00121       FAILED = 2,
00122 // There were no unmasked points
00123       NOGOOD = 3,
00124 // No models set
00125       NOMODELS = 4,
00126 // Number of conditions
00127       nErrorTypes
00128     };
00129 
00130     // Constructor
00131     explicit Fit2D(LogIO& logger);
00132 
00133     // Destructor
00134     ~Fit2D();
00135 
00136     // Copy constructor.  Uses copy semantics except for the logger
00137     // for which a reference copy is made
00138     Fit2D(const Fit2D& other);
00139 
00140     // Assignment operator. Uses copy semantics except for the logger
00141     // for which a reference copy is made
00142     Fit2D& operator=(const Fit2D& other);
00143 
00144     // Add a model to the list to be simultaneously fit and 
00145     // return its index.  Specify the initial guesses for
00146     // the model and a mask indicating whether the parameter
00147     // is fixed (False) during the fit or not.  Returns the
00148     // the model number added (0, 1, 2 etc)
00149     //<group>
00150     uInt addModel (Fit2D::Types type,
00151                    const Vector<Double>& parameters,
00152                    const Vector<Bool>& parameterMask);
00153     uInt addModel(Fit2D::Types type,
00154                    const Vector<Double>& parameters);
00155     //</group>
00156 
00157     // Convert mask from a string to a vector.  The string gives the parameters
00158     // to keep fixed in the fit (f (flux), x (x position), y (y position),
00159     // a (FWHM major axis), b (FWHM minor axis), p (position angle)
00160     static Vector<Bool> convertMask (const String fixedmask,
00161                                      Fit2D::Types type);
00162 
00163 
00164     // Set a pixel selection range.  When the fit is done, only
00165     // pixels in the specified range are included/excluded.
00166     // Only the last call of either of these will be active.
00167     //<group>
00168     void setIncludeRange (Double minVal, Double maxVal);
00169     void setExcludeRange (Double minVal, Double maxVal);
00170     void resetRange();
00171     //</group>
00172 
00173     // Return number of parameters for this type of model
00174     static uInt nParameters (Fit2D::Types type);
00175 
00176     // Recover number of models
00177     uInt nModels() const;
00178 
00179     // Determine an initial estimate for the solution of the specified
00180     // model type to the given data - no compound models are allowable
00181     // in this function.   If you have specified an include
00182     // or exclude pixel range to the fitter, that will be honoured.
00183     // This function does not interact with the addModel function.
00184     // Returns a zero length vector if it fails to make an estimate.
00185     //<group>
00186     Vector<Double> estimate(Fit2D::Types type,
00187                             const MaskedLattice<Float>& data);
00188     Vector<Double> estimate(Fit2D::Types type, const Lattice<Float>& data);
00189     Vector<Double> estimate(Fit2D::Types type, const Array<Float>& data);
00190     Vector<Double> estimate(Fit2D::Types type, const Array<Float>& data,
00191                             const Array<Bool>& mask);
00192     //</group>
00193 
00194     // Do the fit.  Returns an enum value to tell you what happened if the fit failed
00195     // for some reasons.  A message can also be found with function errorMessage if 
00196     // the fit was not successful.  For Array(i,j) i is x and j is y
00197     //<group>
00198     Fit2D::ErrorTypes fit(const MaskedLattice<Float>& data, 
00199                           const Lattice<Float>& sigma);
00200     Fit2D::ErrorTypes fit(const Lattice<Float>& data, 
00201                           const Lattice<Float>& sigma);
00202     Fit2D::ErrorTypes fit(const Array<Float>& data, 
00203                           const Array<Float>& sigma);
00204     Fit2D::ErrorTypes fit(const Array<Float>& data,
00205                           const Array<Bool>& mask, 
00206                           const Array<Float>& sigma);
00207     //</group>
00208 
00209     // Find the residuals to the fit. xOffset and yOffset allow one to provide a data
00210     // array that is offset in space from the grid that was fit. In this way, one
00211     // can fill out a larger image than the subimage that was fit, for example. A negative
00212     // value of xOffset means the supplied data array represents a grid that has a y axis left
00213     // of the grid of pixels that was fit. A negative yOffset value means the supplied data
00214     // array represents a grid that has an x axis that is below the x axis of the grid of pixels
00215     // that was fit.
00216     //<group>
00217     Fit2D::ErrorTypes residual(
00218         Array<Float>& resid, Array<Float>& model,
00219         const Array<Float>& data, Int xOffset=0, int yOffset=0
00220     ) const;
00221     Fit2D::ErrorTypes residual(Array<Float>& resid, Array<Float>& model,
00222                                const MaskedLattice<Float>& data);
00223     Fit2D::ErrorTypes residual(Array<Float>& resid, Array<Float>& model,
00224                                const Lattice<Float>& data);
00225     //</group>
00226     // If function fit failed, you will find a message here
00227     // saying why it failed
00228     String errorMessage () const;
00229 
00230     // Recover solution for either all model components or
00231     // a specific one.  These functions will return an empty vector
00232     // if there is no valid solution.    All available parameters (fixed and
00233     // adjustable) are included in the solution vectors.  
00234     //<group>
00235     Vector<Double> availableSolution () const;
00236     Vector<Double> availableSolution (uInt which) const;
00237     //</group>
00238 
00239     // The errors. All available parameters (fixed and adjustable) are 
00240     // included in the error vectors.  Unsolved for parameters will 
00241     // have error 0.
00242     //<group>
00243     Vector<Double> availableErrors() const;
00244     Vector<Double> availableErrors(uInt which) const;
00245     //</group>
00246 
00247     // The number of iterations that the fitter finished with
00248     uInt numberIterations() const;
00249 
00250     // The chi squared of the fit.  Returns 0 if fit has been done.
00251     Double chiSquared () const;
00252 
00253     // The number of points used for the last fit
00254     uInt numberPoints () const;
00255 
00256     // Return type as a string
00257     static String type(Fit2D::Types type);
00258 
00259     // Return string type as enum (min match)
00260     static Fit2D::Types type(const String& type);
00261 
00262     // Find type of specific model
00263     Fit2D::Types type(uInt which);
00264 
00265     // Convert p.a. (radians) from positive +x -> +y 
00266     // (Fit2D) to positive +y -> -x (Gaussian2D)
00267     static Double paToGauss2D (Double pa) {return pa - C::pi_2;};
00268 
00269     // Convert p.a. (radians) from positive +y -> -x
00270     // (Gaussian2D) to positive +x -> +y (Fit2D)
00271     static Double paFromGauss2D (Double pa) {return pa + C::pi_2;};
00272 
00273 private:
00274 
00275    mutable LogIO itsLogger;
00276    Bool itsValid, itsValidSolution, itsHasSigma;
00277    Bool itsInclude;
00278    Vector<Float> itsPixelRange;
00279    CompoundFunction<AutoDiff<Double> > itsFunction;
00280    NonLinearFitLM<Double> itsFitter;
00281    Vector<Double> itsSolution;
00282    Vector<Double> itsErrors;
00283    Double itsChiSquared;
00284    String itsErrorMessage;
00285    uInt itsNumberPoints;
00286 //
00287    Vector<uInt> itsTypeList;
00288 //
00289    Fit2D::ErrorTypes fitData(const Vector<Double>& values,
00290                              const Matrix<Double>& pos,
00291                              const Vector<Double>& sigma);
00292 
00293 // Returns available (adjustable + fixed) solution for model of
00294 // interest and tells you where it began in the full solution vector
00295 // Does no axial ratio nor position angle conversions from direct
00296 // fit solution vector
00297 // <group>
00298    Vector<Double> availableSolution (uInt& iStart, uInt which) const;
00299    Vector<Double> availableErrors (uInt& iStart, uInt which) const;
00300 // </group>
00301 
00302    Vector<Double> getParams(uInt which) const;
00303    void setParams(const Vector<Double>& params, uInt which);
00304 
00305    Bool includeIt (Float value, const Vector<Float>& range, 
00306                    Int includeIt) const;
00307 
00308    Bool selectData (Matrix<Double>& pos, Vector<Double>& values,
00309                     Vector<Double>& weights,  const Array<Float>& pixels,
00310                     const Array<Bool>& mask, const Array<Float>& sigma);
00311    void piRange (Double& pa) const;
00312 
00313 };
00314 
00315 inline Bool Fit2D::includeIt (Float value, const Vector<Float>& range, 
00316                               Int includeIt) const
00317 {
00318    if (includeIt==0) return True;
00319 //
00320    if (includeIt==1) {
00321       if (value >= range(0) && value <= range(1)) return True;
00322    } else if (value < range(0) || value > range(1)) {
00323       return True;
00324    }
00325 //
00326    return False;
00327 }
00328 
00329 
00330 
00331 
00332 } //# NAMESPACE CASACORE - END
00333 
00334 #endif
00335 
00336 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1