casacore
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 
38 namespace 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 
182 template<class T> class GaussianNDParam : public Function<T>
183 {
184 public:
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
220  virtual ~GaussianNDParam();
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>
247  Vector<T> mean() const;
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 
264 protected:
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.
281 protected:
282  using Function<T>::param_p;
283 public:
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
void setMean(const Vector< T > &mean)
GaussianNDParam(uInt ndim, const T &height)
Vector< T > mean() const
The center ordinate of the Gaussian.
void setVariance(const Vector< T > &variance)
void setCovariance(const Matrix< T > &covar)
void setFlux(const T &flux)
virtual const String & name() const
Give name of function.
GaussianNDParam(uInt ndim, const T &height, const Vector< T > &mean, const Vector< T > &variance)
Matrix< T > covariance() const
The covariance Matrix defines the correlations between all the axes.
void repack(Matrix< T > &covar) const
Functions to convert between internal Vector of parameters and the Covariance Matrix.
void unpack(const Matrix< T > &covar)
GaussianNDParam(const GaussianNDParam< W > &other)
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.
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)
Vector< T > variance() const
The FWHM of the Gaussian is sqrt(8*variance*log(2)).
GaussianNDParam()
Constructs a Gaussian using the indicated height, mean, variance & covariance.
void setHeight(const T &height)
uInt itsDim
dimensionality
virtual ~GaussianNDParam()
Destructor.
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