1 #include <iostream>
2 #include <algorithm>
3 #include "boxm2_ocl_hierarchical_points_to_volume_reg.h"
4 #include <boxm2/ocl/boxm2_ocl_util.h>
5 #include <boct/boct_bit_tree.h>
6 #include <vcl_where_root_dir.h>
7 #include "vnl/vnl_random.h"
8 #ifdef _MSC_VER
9 # include "vcl_msvc_warnings.h"
10 #endif
boxm2_ocl_hierarchical_points_to_volume_reg(boxm2_opencl_cache_sptr & cache,float * pts,boxm2_scene_sptr sceneB,int npts,const bocl_device_sptr & device,bool do_vary_scale)11 boxm2_ocl_hierarchical_points_to_volume_reg::boxm2_ocl_hierarchical_points_to_volume_reg( boxm2_opencl_cache_sptr & cache,
12 float *pts,
13 boxm2_scene_sptr sceneB,
14 int npts,
15 const bocl_device_sptr& device, bool do_vary_scale): boxm2_ocl_reg_points_to_volume_mutual_info(cache,pts,sceneB,device,npts, do_vary_scale), do_vary_scale_(do_vary_scale)
16 {
17
18 }
19
init(vnl_vector<double> const & mu,vnl_vector<double> const & cov)20 bool boxm2_ocl_hierarchical_points_to_volume_reg::init(vnl_vector<double> const& mu, vnl_vector<double> const & cov)
21 {
22 mu_ = mu ;
23 cov_= cov;
24 mu_cost_ = this->cost(mu_);
25 std::cout<<"Mutual Information for current position is "<<mu_cost_<<std::endl;
26 return true;
27 }
28
exhaustive()29 bool boxm2_ocl_hierarchical_points_to_volume_reg::exhaustive()
30 {
31 samples_.clear();
32 int nlevels = 4;
33 unsigned int numsamples_search_width[4] = {3,3,3,3}; // should be odd
34 unsigned int numbestparticales[4] = {1,1,1,1};
35 int params_to_vary = mu_.size() ;
36 if( !do_vary_scale_ )
37 params_to_vary -=1;
38 std::map<double,vnl_vector<double> > samples_sorted;
39 vnl_vector<double> cov = cov_;
40 samples_sorted[mu_cost_] =mu_;
41 for (size_t level = 0; level <nlevels; level++)
42 {
43 std::cout<<"Level #"<<level<<std::endl;
44 mis.clear();
45 samples_.clear();
46 auto iter = samples_sorted.rbegin();
47 for(unsigned int j = 0 ; j < numbestparticales[level] && iter!=samples_sorted.rend(); j++,iter++)
48 {
49 std::cout<<".";
50 std::cout.flush();
51 int searchwidth = numsamples_search_width[level];
52 vnl_vector<double> best_sample = iter->second;
53 int numsamples_per_best_particle = (int) std::pow((float)searchwidth,(float)params_to_vary);
54 for(size_t sampleno = 0 ; sampleno < numsamples_per_best_particle; sampleno++)
55 {
56 int cont=sampleno;
57 vnl_vector<double> curr_sample = best_sample;
58 for(size_t k = 0; k <params_to_vary; k++)
59 {
60 int t = params_to_vary - k - 1;
61 unsigned int var = cont/(unsigned int)std::pow((float)searchwidth,(float)t); // quotient
62 cont = cont - (unsigned int)std::pow((float)searchwidth,(float)t)*var; // remainder
63 double offset = (2*((double)var/((double)searchwidth-1))-1)*cov[t];
64 curr_sample[t] += offset;
65 }
66 double minfo = this->cost(curr_sample, level) ;
67 mis.push_back(minfo);
68 samples_.push_back(curr_sample);
69 }
70 }
71 std::cout<<std::endl;
72 samples_sorted.clear();
73 for(unsigned i = 0 ; i <mis.size(); i++)
74 samples_sorted[mis[i]]=samples_[i];
75 std::cout<<"Mutual Infor for Max Sample"<<this->max_sample()-mu_<<" is "<<this->cost(this->max_sample(),level)<<std::endl;
76 cov = cov_/std::pow((float)2,(float)level+1);
77 }
78 return true;
79 }
max_sample()80 vnl_vector<double> boxm2_ocl_hierarchical_points_to_volume_reg::max_sample()
81 {
82 vnl_vector<double> max_sample( mu_.size(),0.0);
83 double max_pdf = 0.0;
84 for(unsigned i = 0 ; i < mis.size(); i++)
85 if(mis[i] > max_pdf)
86 {
87 max_pdf = mis[i];
88 max_sample = samples_[i];
89 }
90 return max_sample;
91 }
92