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 
35 namespace 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 
178 template<class T> class GaussianND : public GaussianNDParam<T>
179 {
180 public:
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) {}
193  GaussianND(uInt ndim, const T &height) :
194  GaussianNDParam<T>(ndim, height) {}
195  GaussianND(uInt ndim, const T &height, const Vector<T> &mean) :
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.
231 protected:
234 public:
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
const ArgType * FunctionArg
Definition: Function.h:207
Vector< T > mean() const
The center ordinate of the Gaussian.
T height() const
Get or set the peak height of the Gaussian.
GaussianNDParam< T > & operator=(const GaussianNDParam< T > &other)
Copy assignment (deep copy)
virtual uInt ndim() const
Variable dimensionality.
Vector< T > variance() const
The FWHM of the Gaussian is sqrt(8*variance*log(2)).
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