1 // This is core/vpdl/vpdt/vpdt_probability.h
2 #ifndef vpdt_probability_h_
3 #define vpdt_probability_h_
4 //:
5 // \file
6 // \author Matthew Leotta
7 // \brief The basic functions for probability calculations
8 // \date March 13, 2009
9 //
10 // These functions provide default implementations for various probability
11 // calculation functions. They are written in terms of distribution
12 // member functions
13 //
14 // \verbatim
15 // Modifications
16 // <None yet>
17 // \endverbatim
18
19
20 #include <vpdl/vpdt/vpdt_dist_traits.h>
21 #include <vnl/vnl_math.h>
22 #include <cassert>
23 #ifdef _MSC_VER
24 # include <vcl_msvc_warnings.h>
25 #endif
26
27 //: Compute the probability from the density and normalization constant.
28 template <class dist>
29 inline typename vpdt_dist_traits<dist>::scalar_type
vpdt_prob_density(const dist & d,const typename vpdt_dist_traits<dist>::field_type & pt)30 vpdt_prob_density(const dist& d,
31 const typename vpdt_dist_traits<dist>::field_type& pt)
32 {
33 typedef typename vpdt_dist_traits<dist>::scalar_type T;
34 T norm = d.norm_const();
35 if (vnl_math::isinf(norm))
36 return T(0);
37 return norm * d.density(pt);
38 }
39
40
41 //: The probability of being in an axis-aligned box.
42 // The box is defined by two points, the minimum and maximum.
43 // Implemented in terms of \c vpdt_cumulative_prob() by default.
44 template <class dist>
45 typename vpdt_dist_traits<dist>::scalar_type
vpdt_box_prob(const dist & d,const typename vpdt_dist_traits<dist>::field_type & min_pt,const typename vpdt_dist_traits<dist>::field_type & max_pt)46 vpdt_box_prob(const dist& d,
47 const typename vpdt_dist_traits<dist>::field_type& min_pt,
48 const typename vpdt_dist_traits<dist>::field_type& max_pt)
49 {
50 typedef typename vpdt_dist_traits<dist>::scalar_type T;
51 typedef typename vpdt_dist_traits<dist>::field_type F;
52 const unsigned int dim = d.dimension();
53
54 // return zero for ill-defined box
55 for (unsigned int j=0; j<dim; ++j) {
56 if (vpdt_index(max_pt,j)<=vpdt_index(min_pt,j))
57 return T(0);
58 }
59
60 // this method is not tractable for large dimensions
61 assert(sizeof(unsigned long)*8 > dim);
62 // compute the number of corners of the box (2^dim)
63 const unsigned long num_corners = 1 << dim;
64
65 T prob = T(0);
66 F corner(min_pt);
67 for (unsigned long i=0; i<num_corners; ++i)
68 {
69 // In odd dimensions, corners with an odd number of maximal axis
70 // are added to the sum. In even dimensions, corners with an even
71 // number of maximal axis are added. The other corners are subtracted.
72 bool plus = (dim%2 != 1);
73 // create the corner position by selecting elements from max_pt and min_pt
74 for (unsigned int j=0; j<dim; ++j) {
75 bool is_max = (i>>j) & 1;
76 plus ^= is_max; // toggle plus if is_max
77 vpdt_index(corner,j) = is_max?vpdt_index(max_pt,j):vpdt_index(min_pt,j);
78 }
79 if (plus)
80 prob += d.cumulative_prob(corner);
81 else
82 prob -= d.cumulative_prob(corner);
83 }
84
85 return prob;
86 }
87
88 #endif // vpdt_probability_h_
89