1 // This is core/vpdl/vpdl_gaussian_sphere.h
2 #ifndef vpdl_gaussian_sphere_h_
3 #define vpdl_gaussian_sphere_h_
4 //:
5 // \file
6 // \author Matthew Leotta
7 // \date February 11, 2009
8 // \brief A Gaussian with (hyper-)spherical covariance
9 //
10 // \verbatim
11 //  Modifications
12 //   <None yet>
13 // \endverbatim
14 
15 #include <limits>
16 #include "vpdl_gaussian_base.h"
17 #include <vnl/vnl_erf.h>
18 #ifdef _MSC_VER
19 #  include <vcl_msvc_warnings.h>
20 #endif
21 #include <cassert>
22 
23 #include <vpdl/vpdt/vpdt_gaussian.h>
24 #include <vpdl/vpdt/vpdt_probability.h>
25 #include <vpdl/vpdt/vpdt_log_probability.h>
26 
27 //: A Gaussian with (hyper-)spherical covariance
28 template<class T, unsigned int n=0>
29 class vpdl_gaussian_sphere : public vpdl_gaussian_base<T,n>
30 {
31  public:
32   //: the data type used for vectors
33   typedef typename vpdt_field_default<T,n>::type vector;
34   //: the data type used for matrices
35   typedef typename vpdt_field_traits<vector>::matrix_type matrix;
36   //: the type used internally for covariance
37   typedef T covar_type;
38 
39   //: Constructor
40   // Optionally initialize the dimension for when n==0.
41   // Otherwise var_dim is ignored
42   vpdl_gaussian_sphere(unsigned int var_dim = n)
impl_(var_dim)43   : impl_(var_dim) {}
44 
45   //: Constructor - from mean and variance
vpdl_gaussian_sphere(const vector & mean_val,const covar_type & var)46   vpdl_gaussian_sphere(const vector& mean_val, const covar_type& var)
47   : impl_(mean_val,var) {}
48 
49   //: Destructor
50   ~vpdl_gaussian_sphere() override = default;
51 
52   //: Create a copy on the heap and return base class pointer
clone()53   vpdl_distribution<T, n> *clone() const override {
54     return new vpdl_gaussian_sphere<T,n>(*this);
55   }
56 
57   //: Return the run time dimension, which does not equal \c n when \c n==0
dimension()58   unsigned int dimension() const override { return impl_.dimension(); }
59 
60   //: Evaluate the unnormalized density at a point
density(const vector & pt)61   T density(const vector &pt) const override { return impl_.density(pt); }
62 
63   //: Evaluate the probability density at a point
prob_density(const vector & pt)64   T prob_density(const vector &pt) const override {
65     return vpdt_prob_density(impl_,pt);
66   }
67 
68   //: Evaluate the log probability density at a point
log_prob_density(const vector & pt)69   T log_prob_density(const vector &pt) const override {
70     return vpdt_log_prob_density(impl_,pt);
71   };
72 
73   //: Compute the gradient of the unnormalized density at a point
74   // \return the density at \a pt since it is usually needed as well, and
75   //         is often trivial to compute while computing gradient
76   // \retval g the gradient vector
gradient_density(const vector & pt,vector & g)77   T gradient_density(const vector &pt, vector &g) const override {
78     return impl_.gradient_density(pt,g);
79   }
80 
81   //: The normalization constant for the density
82   // When density() is multiplied by this value it becomes prob_density
83   // norm_const() is reciprocal of the integral of density over the entire field
norm_const()84   T norm_const() const override { return impl_.norm_const(); }
85 
86   //: The squared Mahalanobis distance to this point
87   // Non-virtual for efficiency
sqr_mahal_dist(const vector & pt)88   T sqr_mahal_dist(const vector& pt) const
89   {
90     return impl_.sqr_mahal_dist(pt);
91   }
92 
93   //: Evaluate the cumulative distribution function at a point
94   // This is the integral of the density function from negative infinity
95   // (in all dimensions) to the point in question
cumulative_prob(const vector & pt)96   T cumulative_prob(const vector &pt) const override {
97     return impl_.cumulative_prob(pt);
98   }
99 
100   //: The probability of being in an axis-aligned box
101   // The box is defined by two points, the minimum and maximum.
102   // Reimplemented for efficiency since the axis are independent
box_prob(const vector & min_pt,const vector & max_pt)103   T box_prob(const vector &min_pt, const vector &max_pt) const override {
104     const unsigned int dim = this->dimension();
105 
106     double s2 = 1/std::sqrt(2*impl_.covar);
107     // return zero for ill-defined box
108     double prob = T(1);
109     for (unsigned int i=0; i<dim; ++i) {
110       if (vpdt_index(max_pt,i)<=vpdt_index(min_pt,i))
111         return T(0);
112       prob *= (vnl_erf(s2*(vpdt_index(max_pt,i)-vpdt_index(impl_.mean,i))) -
113                vnl_erf(s2*(vpdt_index(min_pt,i)-vpdt_index(impl_.mean,i))))/2;
114     }
115     return static_cast<T>(prob);
116   }
117 
118   //: Access the mean directly
mean()119   const vector &mean() const override { return impl_.mean; }
120 
121   //: Set the mean
set_mean(const vector & mean_val)122   void set_mean(const vector &mean_val) override { impl_.mean = mean_val; }
123 
124   //: Compute the mean of the distribution.
compute_mean(vector & mean_val)125   void compute_mean(vector &mean_val) const override { mean_val = impl_.mean; }
126 
127   //: Access the scalar variance
covariance()128   const covar_type& covariance() const { return impl_.covar; }
129 
130   //: Set the scalar variance
set_covariance(const covar_type & var)131   void set_covariance(const covar_type& var) { impl_.covar = var; }
132 
133   //: Compute the covariance of the distribution.
134   // Should be the identity matrix times var_
compute_covar(matrix & covar)135   void compute_covar(matrix &covar) const override {
136     impl_.compute_covar(covar);
137   }
138 
139  protected:
140   //: the Gaussian implementation from vpdt
141   vpdt_gaussian<vector,covar_type> impl_;
142 };
143 
144 
145 #endif // vpdl_gaussian_sphere_h_
146