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