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