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