casacore
Loading...
Searching...
No Matches
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
39namespace casacore { //# NAMESPACE CASACORE - BEGIN
40
41template<class T> class Lattice;
42template<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
99class Fit2D
100{
101public:
102
103 // Enum describing the different models you can fit
111
112 // Enum describing output error conditions
114// ok
115 OK = 0,
116// Did not converge
118// Solution failed
120// There were no unmasked points
122// No models set
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);
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.
258
259 // The number of points used for the last fit
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
279private:
280
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>
305 Vector<Double> availableErrors (uInt& iStart, uInt which) const;
306// </group>
307
309 void setParams(const Vector<Double>& params, uInt which);
310
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
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 Array< T > &data, const Array< Bool > &mask, const Array< T > &sigma)
Fit2D::ErrorTypes residual(Array< Float > &resid, Array< Float > &model, const Lattice< Float > &data)
Vector< Double > availableErrors(uInt which) const
uInt itsNumberPoints
Definition Fit2D.h:291
LogIO itsLogger
Definition Fit2D.h:281
Vector< Double > availableSolution() const
Recover solution for either all model components or a specific one.
Bool selectData(Matrix< Double > &pos, Vector< Double > &values, Vector< Double > &weights, const Array< T > &pixels, const Array< Bool > &mask, const Array< T > &sigma)
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.
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
Vector< Double > availableSolution(uInt &iStart, uInt which) const
Returns available (adjustable + fixed) solution for model of interest and tells you where it began in...
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.
Fit2D & operator=(const Fit2D &other)
Assignment operator.
Bool itsValid
Definition Fit2D.h:282
Vector< Double > estimate(Fit2D::Types type, const Array< T > &data)
Vector< uInt > itsTypeList
Definition Fit2D.h:293
Double chiSquared() const
The chi squared of the fit.
Fit2D::ErrorTypes fit(const MaskedLattice< T > &data, const Lattice< T > &sigma)
Do the fit.
Vector< Double > estimate(Fit2D::Types type, const Lattice< T > &data)
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)
Vector< Double > availableErrors(uInt &iStart, uInt which) const
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)
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)
Double itsChiSquared
Definition Fit2D.h:289
Types
Enum describing the different models you can fit.
Definition Fit2D.h:104
Vector< Double > availableErrors() const
The errors.
static Vector< Bool > convertMask(const String fixedmask, Fit2D::Types type)
Convert mask from a string to a vector.
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.
Vector< Double > availableSolution(uInt which) const
Vector< Double > itsErrors
Definition Fit2D.h:288
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 > estimate(Fit2D::Types type, const Array< T > &data, const Array< Bool > &mask)
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.
Vector< Double > getParams(uInt which) const
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