1 #include <iostream>
2 #include <algorithm>
3 #include <cmath>
4 #include "bvpl_gauss3d_kernel_factory.h"
5 //:
6 // \file
7
8 #ifdef _MSC_VER
9 # include "vcl_msvc_warnings.h"
10 #endif
11 #include "vnl/vnl_math.h"
12 #include "vnl/vnl_float_3.h"
13 #include <bsta/bsta_gauss_if3.h>
14
15 // Default Constructor
bvpl_gauss3d_kernel_factory()16 bvpl_gauss3d_kernel_factory::bvpl_gauss3d_kernel_factory()
17 {
18 sigma1_ = 0.0f;
19 sigma2_ = 0.0f;
20 sigma3_ = 0.0f;
21 angular_resolution_ = 0;
22 rotation_axis_ = canonical_rotation_axis_;
23 angle_ = 0.0f;
24 }
25
26 //: Constructs a kernel form gaussian spheroid with sigma parameters s1 and s2. i.e. Cov is diagonal with entries s1, s2, s2
bvpl_gauss3d_kernel_factory(float s1,float s2)27 bvpl_gauss3d_kernel_factory::bvpl_gauss3d_kernel_factory(float s1, float s2)
28 {
29 //set variances of kernel
30 sigma1_ = s1;
31 sigma2_ = s2;
32 sigma3_ = s2;
33
34 //this skernel is symmetric around main axis
35 angular_resolution_=0;
36
37 //initialize variables
38 angle_ = 0.0f;
39 rotation_axis_ = canonical_rotation_axis_;
40 parallel_axis_ = canonical_parallel_axis_;
41
42 //create the default kernel
43 create_canonical();
44 }
45
46 //: Constructs a kernel form gaussian ellipsoid with sigma parameters s1, s2 and s3. i.e. Cov is diagonal with entries s1, s2, s3
bvpl_gauss3d_kernel_factory(float s1,float s2,float s3,float supp1,float supp2,float supp3)47 bvpl_gauss3d_kernel_factory::bvpl_gauss3d_kernel_factory(float s1, float s2, float s3,float supp1 , float supp2 , float supp3 )
48 {
49 //set variances of kernel
50 sigma1_ = s1;
51 sigma2_ = s2;
52 sigma3_ = s3;
53 supp1_ = supp1;
54 supp2_ = supp2;
55 supp3_ = supp3;
56
57 //this value is a meant as a limit there is not theoretical meaning to it
58 angular_resolution_= float(vnl_math::pi/16.0);
59
60 //initialize variables
61 angle_ = 0.0f;
62 rotation_axis_ = canonical_rotation_axis_;
63 parallel_axis_ = canonical_parallel_axis_;
64
65
66 //create the default kernel
67 create_canonical();
68 }
69
create_canonical()70 void bvpl_gauss3d_kernel_factory::create_canonical()
71 {
72 bsta_gauss_if3 gauss_kernel(vnl_float_3(0,0,0), vnl_float_3(sigma1_*sigma1_, sigma2_*sigma2_, sigma3_*sigma3_));
73
74 typedef vgl_point_3d<float> point_3d;
75 typedef bvpl_kernel_dispatch dispatch;
76
77 //This is the support of the kernel
78 float min_x = -std::floor(supp1_*sigma1_+0.01f);
79 float max_x = std::floor(supp1_*sigma1_+0.01f);
80 float min_y = -std::floor(supp2_*sigma2_+0.01f);
81 float max_y = std::floor(supp2_*sigma2_+0.01f);
82 float min_z = -std::floor(supp3_*sigma3_+0.01f);
83 float max_z = std::floor(supp3_*sigma3_+0.01f);
84 float l1_norm = 0.0f;
85
86 for (float x=min_x; x<= max_x; x+=1.f)
87 {
88 for (float y= min_y; y<= max_y; y+=1.f)
89 {
90 for (float z= min_z; z<= max_z; z+=1.f)
91 {
92 vnl_float_3 pt(x,y,z);
93 float val = gauss_kernel.prob_density(pt);
94 canonical_kernel_.emplace_back(point_3d(x,y,z), dispatch(val));
95 l1_norm +=val;
96 }
97 }
98 }
99
100 //normalize to L1 norm
101 auto k_it = canonical_kernel_.begin();
102 float norm = 0.0;
103 for(; k_it != canonical_kernel_.end(); k_it++)
104 {
105 k_it->second.c_ = k_it->second.c_ / l1_norm;
106 norm+=k_it->second.c_;
107 }
108 std::cout << "Canonical kernel should have l1_norm, norm= " << norm <<std::endl;
109
110 //set the dimension of the 3-d grid
111 max_point_.set(int(max_x),int(max_y),int(max_z));
112 min_point_.set(int(min_x),int(min_y),int(min_z));
113
114 //set the current kernel
115 kernel_ = canonical_kernel_;
116 factory_name_ = name();
117 }
118