casacore
Loading...
Searching...
No Matches
GaussianNDParam.h
Go to the documentation of this file.
1//# GaussianNDParam.h: Multidimensional Gaussian class parameters
2//# Copyright (C) 2001,2002,2004,2005
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 SCIMATH_GAUSSIANNDPARAM_H
29#define SCIMATH_GAUSSIANNDPARAM_H
30
31//# Includes
32#include <casacore/casa/aips.h>
33#include <casacore/scimath/Functionals/Function.h>
34#include <casacore/casa/Arrays/Matrix.h>
35#include <casacore/casa/Arrays/Vector.h>
36#include <casacore/casa/BasicSL/String.h>
37
38namespace casacore { //# NAMESPACE CASACORE - BEGIN
39
40// <summary> A Multi-dimensional Gaussian parameter handling. </summary>
41
42// <use visibility=local>
43
44// <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tGaussianND" demos="dGaussianND">
45// </reviewed>
46
47// <prerequisite>
48// <li> <linkto class="FunctionParam">FunctionParam</linkto> class
49// <li> <linkto class="Function">Function</linkto> class
50// </prerequisite>
51
52// <synopsis>
53// A <src>GaussianND</src> is used to calculate Gaussian functions of any
54// dimension. A <linkto class=Gaussian1D> Gaussian1D </linkto> class exists
55// which is more appropriate for one dimensional Gaussian functions, and a
56// <linkto class=Gaussian2D> Gaussian2D </linkto> class exists for two
57// dimensional functions.
58//
59// A statistical description of the multi-dimensional Gaussian is used (see
60// Kendall & Stuart "The Advanced Theory of Statistics"). A Gaussian is
61// defined in terms of its height, mean (which is the location of the peak
62// value), variance, (a measure of the width of the Gaussian), and
63// covariance which skews the distribution with respect to the Axes.
64//
65// In the general description the variance and covariance are specified
66// using a covariance matrix. This is defined as (for a 4 dimensional
67// Gaussian):
68// <srcblock>
69// V = | s1*s1 r12*s1*s2 r13*s1*s3 r14*s1*s4 |
70// | r12*s1*s2 s2*s2 r23*s2*s3 r24*s2*s4 |
71// | r13*s1*s3 r23*s2*s3 s3*s3 r34*s3*s4 |
72// | r14*s1*s4 r24*s2*s4 r34*s3*s4 s4*s4 |
73// </srcblock>
74// where s1 (<src>sigma1</src>) is the standard deviation of the Gaussian with
75// respect to the first axis, and r12 (<src>rho12</src>) is the correlation
76// between the the first and second axis. The correlation MUST be between -1
77// and 1, and this class checks this as well as ensuring that the diagonal
78// is positive.
79//
80// <note role=warning> It is possible to have symmetric matrices that are of
81// the above described form (ie. symmetric with <src>-1 <= rho(ij) <=1</src>)
82// that do
83// not generate a Gaussian function. This is because the Matrix is NOT
84// positive definite (The limits on <src>rho(ij)</src> are upper limits).
85// This class
86// does check that the covariance Matrix is positive definite and will throw
87// an exception (AipsError) if it is not.</note>
88//
89// The covariance Matrix can be specified by only its upper or lower
90// triangular regions (ie. with zeros in the other triangle), otherwise it
91// MUST be symmetric.
92//
93// The Gaussian that is constructed from this covariance Matrix (V), along
94// with mean (u) and height (h) is:
95// <srcblock>
96// f(x) = h*exp( -1/2 * (x-u) * V^(-1) * (x-u))
97// </srcblock>
98// where x, and u are vectors whose length is the dimensionality of the
99// Gaussian and V^(-1) is the inverse of the covariance Matrix defined
100// above. For a two dimensional Gaussian with zero mean this expression
101// reduces to:
102// <srcblock>
103// f(x) = h*exp(-1/(2*(1-r12^2))*(x1^2/s1^2 - 2*r12*x1*x2/(s1*s2) + x2^2/s2^2))
104// </srcblock>
105//
106// The amplitude of the Gaussian can be defined in two ways, either using
107// the peak height (as is done in the constructors, and the setHeight
108// function) or using the setFlux function. The flux in this context is the
109// analytic integral of the Gaussian over all dimensions. Using the setFlux
110// function does not modify the shape of the Gaussian just its height.
111//
112// All the parameters of the Gaussian except its dimensionality can be
113// modified using the set/get functions.
114//
115// The parameter interface (see
116// <linkto class="FunctionParam">FunctionParam</linkto> class),
117// is used to provide an interface to the
118// <linkto module="Fitting"> Fitting </linkto> classes.
119// There are always 4 parameter sets.
120// <note role=warning> Note that the actual variance/covariance
121// parameters are the inverse matrix of the variance/covariance matrix given
122// by the user</note>.
123// The actual parameters are in order:
124// <ol>
125// <li> height (1 term). No assumptions on what quantity the height
126// represents, and it can be negative (enumerated by <src>HEIGHT</src>)
127// <li> mean (ndim terms) (enumerated by <src>CENTER</src>).
128// <li> variance (ndim terms). The variance is always positive, and an
129// exception (AipsError) will be thrown if you try to set a negative
130// value.
131// <li> covariance (ndim*(ndim-1)/2 terms) The order is (assuming ndim=5)
132// v12,v13,v14,v15,v23,v24,v25,v34,v35,v45. The restrictions described
133// above for the covariance (ie. -1 < r12 < +1) are enforced.
134// </ol>
135// </synopsis>
136
137// <example>
138// Construct a two dimensional Gaussian with mean=(0,1), variance=(.1,7) and
139// height = 1;
140// <srcblock>
141// uInt ndim = 2;
142// Float height = 1;
143// Vector<Float> mean(ndim); mean(0) = 0, mean(1) = 1;
144// Vector<Float> variance(ndim); variance(0) = .1, variance(1) = 7;
145// GaussianND<Float> g(ndim, height, mean, variance);
146// Vector<Float> x(ndim); x = 0;
147// cout << "g("<< x <<") = " << g(x) <<endl; // g([0,0])=1*exp(-1/2*1/7);
148// x(1)++;
149// cout << "g("<< x <<") = " <<g(x) <<endl; // g([0,1])= 1
150// cout << "Height: " << g.height() <<endl; // Height: 1
151// cout << "Flux: " << g.flux() << endl; // Flux: 2*Pi*Sqrt(.1*7)
152// cout << "Mean: " << g.mean() << endl; // Mean: [0, -1]
153// cout << "Variance: " << g.variance() <<endl; // Variance: [.1, 7]
154// cout << "Covariance: "<< g.covariance()<<endl;// Covariance: [.1, 0]
155// // [0, 7]
156// g.setFlux(1);
157// cout << "g("<< x <<") = " <<g(x) <<endl; //g([0,1])=1/(2*Pi*Sqrt(.7))
158// cout << "Height: " << g.height() <<endl; // Height: 1/(2*Pi*Sqrt(.7))
159// cout << "Flux: " << g.flux() << endl; // Flux: 1
160// cout << "Mean: " << g.mean() << endl; // Mean: [0, -1]
161// cout << "Variance: " << g.variance() <<endl; // Variance: [.1, 7]
162// cout << "Covariance: "<< g.covariance()<<endl;// Covariance: [.1, 0]
163// // [0, 7]
164// </srcblock>
165// </example>
166
167// <motivation>
168// A Gaussian Functional was needed for modeling the sky with a series of
169// components. It was later realised that it was too general and Gaussian2D
170// was written.
171// </motivation>
172
173// <templating arg=T>
174// <li> T should have standard numerical operators and exp() function. Current
175// implementation only tested for real types.
176// </templating>
177
178// <todo asof="2001/08/19">
179// <li> Nothing I know off, apart from possible optimization
180// </todo>
181
182template<class T> class GaussianNDParam : public Function<T>
183{
184public:
185 //# Enumerations
186 enum { HEIGHT=0, CENTER};
187
188 //# Constructors
189 // Constructs a Gaussian using the indicated height, mean, variance &
190 // covariance.
191 // ndim defaults to 2,
192 // mean defaults to 0,
193 // height to <src> Pi^(-ndim/2)</src> (the flux is unity)
194 // variance defaults to 1.0,
195 // covariance defaults to 0.0,
196 // <group>
202 const Vector<T> &variance);
204 const Matrix<T> &covar);
205 // </group>
206
207 // Copy constructor (deep copy)
208 // <group>
210 template <class W>
212 Function<T>(other),
213 itsDim(other.itsDim), itsFlux2Hgt(other.itsFlux2Hgt) {}
214 // </group>
215
216 // Copy assignment (deep copy)
218
219 // Destructor
221
222 //# Operators
223
224 //# Member functions
225 // Give name of function
226 virtual const String &name() const { static String x("gaussiannd");
227 return x; }
228
229 // Variable dimensionality
230 virtual uInt ndim() const { return itsDim; }
231
232 // Get or set the peak height of the Gaussian
233 // <group>
234 T height() const { return param_p[HEIGHT]; }
235 void setHeight(const T &height) { param_p[HEIGHT] = height; }
236 // </group>
237
238 // The analytical integrated area underneath the Gaussian. Use these
239 // functions as an alternative to the height functions.
240 // <group>
241 T flux() const;
242 void setFlux(const T &flux);
243 // </group>
244
245 // The center ordinate of the Gaussian
246 // <group>
248 void setMean(const Vector<T> &mean);
249 // </group>
250
251 // The FWHM of the Gaussian is <src>sqrt(8*variance*log(2))</src>.
252 // The variance MUST be positive
253 // <group>
256 // </group>
257
258 //The covariance Matrix defines the correlations between all the axes.
259 //<group>
261 void setCovariance(const Matrix<T> &covar);
262 // </group>
263
264protected:
265 //# Data
266 // dimensionality
268 // factor to convert from flux to height
270
271 //# Methods
272 // Functions to convert between internal Vector of parameters
273 // and the Covariance
274 // Matrix
275 // <group>
276 void repack(Matrix<T> &covar) const;
277 void unpack(const Matrix<T> &covar);
278 // </group>
279
280 //# Make members of parent classes known.
281protected:
282 using Function<T>::param_p;
283public:
284 using Function<T>::nparameters;
285};
286
287
288} //# NAMESPACE CASACORE - END
289
290#ifndef CASACORE_NO_AUTO_TEMPLATES
291#include <casacore/scimath/Functionals/GaussianNDParam.tcc>
292#endif //# CASACORE_NO_AUTO_TEMPLATES
293#endif
FunctionParam< T > param_p
The parameters and masks.
Definition Function.h:332
uInt nparameters() const
Returns the number of parameters.
Definition Function.h:230
Vector< T > variance() const
The FWHM of the Gaussian is sqrt(8*variance*log(2)).
void setMean(const Vector< T > &mean)
GaussianNDParam(uInt ndim, const T &height)
Matrix< T > covariance() const
The covariance Matrix defines the correlations between all the axes.
void setVariance(const Vector< T > &variance)
void setCovariance(const Matrix< T > &covar)
void setFlux(const T &flux)
GaussianNDParam(uInt ndim, const T &height, const Vector< T > &mean, const Vector< T > &variance)
void repack(Matrix< T > &covar) const
Functions to convert between internal Vector of parameters and the Covariance Matrix.
virtual const String & name() const
Give name of function.
void unpack(const Matrix< T > &covar)
GaussianNDParam(const GaussianNDParam< W > &other)
T height() const
Get or set the peak height of the Gaussian.
virtual uInt ndim() const
Variable dimensionality.
GaussianNDParam(const GaussianNDParam &other)
Copy constructor (deep copy)
T itsFlux2Hgt
factor to convert from flux to height
T flux() const
The analytical integrated area underneath the Gaussian.
GaussianNDParam(uInt ndim, const T &height, const Vector< T > &mean, const Matrix< T > &covar)
GaussianNDParam(uInt ndim, const T &height, const Vector< T > &mean)
GaussianNDParam< T > & operator=(const GaussianNDParam< T > &other)
Copy assignment (deep copy)
GaussianNDParam()
Constructs a Gaussian using the indicated height, mean, variance & covariance.
void setHeight(const T &height)
uInt itsDim
dimensionality
virtual ~GaussianNDParam()
Destructor.
Vector< T > mean() const
The center ordinate of the Gaussian.
String: the storage and methods of handling collections of characters.
Definition String.h:225
this file contains all the compiler specific defines
Definition mainpage.dox:28
unsigned int uInt
Definition aipstype.h:51