casacore
Fit2D.h
Go to the documentation of this file.
1 //# Fit2D.h: Class to fit 2-D objects to Lattices or Arrays
2 //# Copyright (C) 1997,1998,1999,2000,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 //#
26 //# $Id$
27 
28 #ifndef LATTICES_FIT2D_H
29 #define LATTICES_FIT2D_H
30 
31 //# Includes
32 #include <casacore/casa/aips.h>
33 #include <casacore/casa/Arrays/ArrayFwd.h>
34 #include <casacore/scimath/Functionals/CompoundFunction.h>
35 #include <casacore/casa/BasicSL/Constants.h>
36 #include <casacore/scimath/Fitting/NonLinearFitLM.h>
37 #include <casacore/casa/Logging/LogIO.h>
38 
39 namespace casacore { //# NAMESPACE CASACORE - BEGIN
40 
41 template<class T> class Lattice;
42 template<class T> class MaskedLattice;
43 
44 
45 // <summary>
46 // Fit 2-D objects to 2-D Lattices or Arrays
47 // </summary>
48 
49 // <use visibility=export>
50 
51 // <reviewed reviewer="" date="" tests="">
52 // </reviewed>
53 
54 // <prerequisite>
55 // <li> <linkto class=Lattice>Lattice</linkto>
56 // </prerequisite>
57 
58 // <synopsis>
59 // This class allows you to fit different types of 2-D models
60 // to either Lattices or Arrays. These must be 2 dimensional;
61 // for Lattices, the appropriate 2-D Lattice can be made with
62 // the SubLattice class.
63 //
64 // You may fit more than one model simultaneously to the data.
65 // Models are added with the addModel method. With this method,
66 // you also specify the initial guesses of the parameters of
67 // the model. Any parameters involving coordinates are
68 // expected in zero-relative absolute pixel coordinates (e.g. the centre of
69 // a model). Additionally with the addModel method,
70 // you may specify which parameters are to be held fixed
71 // during the fitting process. This is done with the
72 // parameterMask Vector which is in the same order as the
73 // parameter Vector. A value of True indicates the parameter
74 // will be fitted for. Presently, when you say fix the minor axis,
75 // you really end up fixing the axial ratio (internals). I don't
76 // have a solution for this presently.
77 //
78 // For Gaussians, the parameter Vector (input or output) consists, in order, of
79 // the peak, x location, y location, FWHM of major axis, FWHM of minor axis,
80 // and position angle of the major axis (in radians). The
81 // position angle is positive +x to +y
82 // in the pixel coordinate system ([0,0] in center of image) and
83 // in the range -2pi to 2pi. When the solution is recovered, the
84 // position angle will be in the range 0 to pi.
85 //
86 // </synopsis>
87 // <example>
88 // <srcblock>
89 // </srcblock>
90 // </example>
91 
92 // <todo asof="1998/12/11">
93 // <li> template it
94 // <li> Speed up some Array calculations indexed with IPositions
95 // <li> Don't handle Lattices simply by getting pixels into Arrays
96 // <li> make an addModel interface taking functionals
97 // </todo>
98 
99 class Fit2D
100 {
101 public:
102 
103  // Enum describing the different models you can fit
104  enum Types {
105  GAUSSIAN = 0,
106  DISK = 1,
107  LEVEL=2,
108  PLANE=3,
109  nTypes
110  };
111 
112  // Enum describing output error conditions
113  enum ErrorTypes {
114 // ok
115  OK = 0,
116 // Did not converge
118 // Solution failed
119  FAILED = 2,
120 // There were no unmasked points
121  NOGOOD = 3,
122 // No models set
123  NOMODELS = 4,
124 // Number of conditions
126  };
127 
128  // Constructor
129  explicit Fit2D(LogIO& logger);
130 
131  // Destructor
133 
134  // Copy constructor. Uses copy semantics except for the logger
135  // for which a reference copy is made
136  Fit2D(const Fit2D& other);
137 
138  // Assignment operator. Uses copy semantics except for the logger
139  // for which a reference copy is made
140  Fit2D& operator=(const Fit2D& other);
141 
142  // Add a model to the list to be simultaneously fit and
143  // return its index. Specify the initial guesses for
144  // the model and a mask indicating whether the parameter
145  // is fixed (False) during the fit or not. Returns the
146  // the model number added (0, 1, 2 etc)
147  //<group>
149  const Vector<Double>& parameters,
150  const Vector<Bool>& parameterMask);
152  const Vector<Double>& parameters);
153  //</group>
154 
155  // Convert mask from a string to a vector. The string gives the parameters
156  // to keep fixed in the fit (f (flux), x (x position), y (y position),
157  // a (FWHM major axis), b (FWHM minor axis), p (position angle)
158  static Vector<Bool> convertMask (const String fixedmask,
160 
161 
162  // Set a pixel selection range. When the fit is done, only
163  // pixels in the specified range are included/excluded.
164  // Only the last call of either of these will be active.
165  //<group>
166  void setIncludeRange (Double minVal, Double maxVal);
167  void setExcludeRange (Double minVal, Double maxVal);
168  void resetRange();
169  //</group>
170 
171  // Return number of parameters for this type of model
173 
174  // Recover number of models
175  uInt nModels() const;
176 
177  // Determine an initial estimate for the solution of the specified
178  // model type to the given data - no compound models are allowable
179  // in this function. If you have specified an include
180  // or exclude pixel range to the fitter, that will be honoured.
181  // This function does not interact with the addModel function.
182  // Returns a zero length vector if it fails to make an estimate.
183  //<group>
185  const MaskedLattice<T>& data);
186  template<class T> Vector<Double> estimate(
187  Fit2D::Types type, const Lattice<T>& data
188  );
189  template<class T> Vector<Double> estimate(
190  Fit2D::Types type, const Array<T>& data
191  );
192  template<class T> Vector<Double> estimate(
193  Fit2D::Types type, const Array<T>& data, const Array<Bool>& mask
194  );
195  //</group>
196 
197  // Do the fit. Returns an enum value to tell you what happened if the fit failed
198  // for some reasons. A message can also be found with function errorMessage if
199  // the fit was not successful. For Array(i,j) i is x and j is y
200  //<group>
201  template <class T> Fit2D::ErrorTypes fit(const MaskedLattice<T>& data,
202  const Lattice<T>& sigma);
203  template <class T> Fit2D::ErrorTypes fit(const Lattice<T>& data,
204  const Lattice<T>& sigma);
205  template <class T> Fit2D::ErrorTypes fit(const Array<T>& data,
206  const Array<T>& sigma);
207  template <class T> Fit2D::ErrorTypes fit(const Array<T>& data,
208  const Array<Bool>& mask,
209  const Array<T>& sigma);
210  //</group>
211 
212  // Find the residuals to the fit. xOffset and yOffset allow one to provide a data
213  // array that is offset in space from the grid that was fit. In this way, one
214  // can fill out a larger image than the subimage that was fit, for example. A negative
215  // value of xOffset means the supplied data array represents a grid that has a y axis left
216  // of the grid of pixels that was fit. A negative yOffset value means the supplied data
217  // array represents a grid that has an x axis that is below the x axis of the grid of pixels
218  // that was fit.
219  // NOTE these may need to be templated at some point in the future. My
220  // current need does not require they be templated. - dmehring 29jun2018
221  //<group>
222  template <class T> Fit2D::ErrorTypes residual(
223  Array<T>& resid, Array<T>& model,
224  const Array<T>& data, Int xOffset=0, int yOffset=0
225  ) const;
226 
228  const MaskedLattice<Float>& data);
230  const Lattice<Float>& data);
231  //</group>
232  // If function fit failed, you will find a message here
233  // saying why it failed
235 
236  // Recover solution for either all model components or
237  // a specific one. These functions will return an empty vector
238  // if there is no valid solution. All available parameters (fixed and
239  // adjustable) are included in the solution vectors.
240  //<group>
243  //</group>
244 
245  // The errors. All available parameters (fixed and adjustable) are
246  // included in the error vectors. Unsolved for parameters will
247  // have error 0.
248  //<group>
251  //</group>
252 
253  // The number of iterations that the fitter finished with
255 
256  // The chi squared of the fit. Returns 0 if fit has been done.
257  Double chiSquared () const;
258 
259  // The number of points used for the last fit
260  uInt numberPoints () const;
261 
262  // Return type as a string
264 
265  // Return string type as enum (min match)
266  static Fit2D::Types type(const String& type);
267 
268  // Find type of specific model
270 
271  // Convert p.a. (radians) from positive +x -> +y
272  // (Fit2D) to positive +y -> -x (Gaussian2D)
273  static Double paToGauss2D (Double pa) {return pa - C::pi_2;};
274 
275  // Convert p.a. (radians) from positive +y -> -x
276  // (Gaussian2D) to positive +x -> +y (Fit2D)
277  static Double paFromGauss2D (Double pa) {return pa + C::pi_2;};
278 
279 private:
280 
281  mutable LogIO itsLogger;
292 
294 
296  const Matrix<Double>& pos,
297  const Vector<Double>& sigma);
298 
299 // Returns available (adjustable + fixed) solution for model of
300 // interest and tells you where it began in the full solution vector
301 // Does no axial ratio nor position angle conversions from direct
302 // fit solution vector
303 // <group>
304  Vector<Double> availableSolution (uInt& iStart, uInt which) const;
305  Vector<Double> availableErrors (uInt& iStart, uInt which) const;
306 // </group>
307 
309  void setParams(const Vector<Double>& params, uInt which);
310 
311  Bool includeIt (Double value, const Vector<Double>& range,
312  Int includeIt) const;
313 
314  template <class T> Bool selectData (
315  Matrix<Double>& pos, Vector<Double>& values,
316  Vector<Double>& weights, const Array<T>& pixels,
317  const Array<Bool>& mask, const Array<T>& sigma
318  );
319 
320  void piRange (Double& pa) const;
321 
322 };
323 
325  Int includeIt) const
326 {
327  if (includeIt==0) return True;
328 //
329  if (includeIt==1) {
330  if (value >= range(0) && value <= range(1)) return True;
331  } else if (value < range(0) || value > range(1)) {
332  return True;
333  }
334 //
335  return False;
336 }
337 
338 } //# NAMESPACE CASACORE - END
339 
340 #ifndef CASACORE_NO_AUTO_TEMPLATES
341 #include <casacore/lattices/LatticeMath/Fit2D2.tcc>
342 #endif //# CASACORE_NO_AUTO_TEMPLATES
343 
344 #endif
345 
346 
Fit2D::ErrorTypes fit(const Array< T > &data, const Array< Bool > &mask, const Array< T > &sigma)
Fit2D::ErrorTypes residual(Array< Float > &resid, Array< Float > &model, const Lattice< Float > &data)
uInt itsNumberPoints
Definition: Fit2D.h:291
Vector< Double > estimate(Fit2D::Types type, const Lattice< T > &data)
LogIO itsLogger
Definition: Fit2D.h:277
Bool selectData(Matrix< Double > &pos, Vector< Double > &values, Vector< Double > &weights, const Array< T > &pixels, const Array< Bool > &mask, const Array< T > &sigma)
Vector< Double > estimate(Fit2D::Types type, const Array< T > &data)
Vector< Double > availableSolution(uInt &iStart, uInt which) const
Returns available (adjustable + fixed) solution for model of interest and tells you where it began in...
uInt addModel(Fit2D::Types type, const Vector< Double > &parameters, const Vector< Bool > &parameterMask)
Add a model to the list to be simultaneously fit and return its index.
String errorMessage() const
If function fit failed, you will find a message here saying why it failed.
Vector< Double > estimate(Fit2D::Types type, const MaskedLattice< T > &data)
Determine an initial estimate for the solution of the specified model type to the given data - no com...
Fit2D::ErrorTypes fit(const Lattice< T > &data, const Lattice< T > &sigma)
Fit2D(LogIO &logger)
Constructor.
static Double paToGauss2D(Double pa)
Convert p.a.
Definition: Fit2D.h:273
Fit2D(const Fit2D &other)
Copy constructor.
Vector< Double > itsSolution
Definition: Fit2D.h:287
void setIncludeRange(Double minVal, Double maxVal)
Set a pixel selection range.
NonLinearFitLM< Double > itsFitter
Definition: Fit2D.h:286
static uInt nParameters(Fit2D::Types type)
Return number of parameters for this type of model.
Bool itsValid
Definition: Fit2D.h:282
Vector< Double > availableSolution() const
Recover solution for either all model components or a specific one.
Vector< uInt > itsTypeList
Definition: Fit2D.h:293
Vector< Double > availableErrors() const
The errors.
Double chiSquared() const
The chi squared of the fit.
Fit2D::ErrorTypes fit(const MaskedLattice< T > &data, const Lattice< T > &sigma)
Do the fit.
static Vector< Bool > convertMask(const String fixedmask, Fit2D::Types type)
Convert mask from a string to a vector.
CompoundFunction< AutoDiff< Double > > itsFunction
Definition: Fit2D.h:285
static String type(Fit2D::Types type)
Return type as a string.
void setExcludeRange(Double minVal, Double maxVal)
String itsErrorMessage
Definition: Fit2D.h:290
uInt addModel(Fit2D::Types type, const Vector< Double > &parameters)
Bool itsHasSigma
Definition: Fit2D.h:282
Fit2D::ErrorTypes residual(Array< Float > &resid, Array< Float > &model, const MaskedLattice< Float > &data)
Vector< Double > availableErrors(uInt &iStart, uInt which) const
uInt numberIterations() const
The number of iterations that the fitter finished with.
Fit2D::ErrorTypes fitData(const Vector< Double > &values, const Matrix< Double > &pos, const Vector< Double > &sigma)
Vector< Double > estimate(Fit2D::Types type, const Array< T > &data, const Array< Bool > &mask)
Vector< Double > getParams(uInt which) const
Double itsChiSquared
Definition: Fit2D.h:289
Types
Enum describing the different models you can fit.
Definition: Fit2D.h:104
void piRange(Double &pa) const
static Fit2D::Types type(const String &type)
Return string type as enum (min match)
Fit2D::Types type(uInt which)
Find type of specific model.
uInt nModels() const
Recover number of models.
Fit2D & operator=(const Fit2D &other)
Assignment operator.
Vector< Double > itsErrors
Definition: Fit2D.h:288
Vector< Double > availableErrors(uInt which) const
void setParams(const Vector< Double > &params, uInt which)
Bool itsInclude
Definition: Fit2D.h:283
Fit2D::ErrorTypes fit(const Array< T > &data, const Array< T > &sigma)
ErrorTypes
Enum describing output error conditions.
Definition: Fit2D.h:113
@ FAILED
Solution failed.
Definition: Fit2D.h:119
@ NOMODELS
No models set.
Definition: Fit2D.h:123
@ NOGOOD
There were no unmasked points.
Definition: Fit2D.h:121
@ NOCONVERGE
Did not converge.
Definition: Fit2D.h:117
@ nErrorTypes
Number of conditions.
Definition: Fit2D.h:125
Bool itsValidSolution
Definition: Fit2D.h:282
Vector< Double > availableSolution(uInt which) const
Vector< Double > itsPixelRange
Definition: Fit2D.h:284
Fit2D::ErrorTypes residual(Array< T > &resid, Array< T > &model, const Array< T > &data, Int xOffset=0, int yOffset=0) const
Find the residuals to the fit.
uInt numberPoints() const
The number of points used for the last fit.
Bool includeIt(Double value, const Vector< Double > &range, Int includeIt) const
Definition: Fit2D.h:324
~Fit2D()
Destructor.
static Double paFromGauss2D(Double pa)
Convert p.a.
Definition: Fit2D.h:277
String: the storage and methods of handling collections of characters.
Definition: String.h:225
const Double pi_2
pi/2
this file contains all the compiler specific defines
Definition: mainpage.dox:28
LatticeExprNode pa(const LatticeExprNode &left, const LatticeExprNode &right)
This function finds 180/pi*atan2(left,right)/2.
const Bool False
Definition: aipstype.h:44
unsigned int uInt
Definition: aipstype.h:51
LatticeExprNode mask(const LatticeExprNode &expr)
This function returns the mask of the given expression.
int Int
Definition: aipstype.h:50
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
LatticeExprNode value(const LatticeExprNode &expr)
This function returns the value of the expression without a mask.
const Bool True
Definition: aipstype.h:43
double Double
Definition: aipstype.h:55