1 #include <iostream>
2 #include <algorithm>
3 #include "boxm2_ocl_hierarchical_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_reg(boxm2_opencl_cache_sptr & cache,boxm2_scene_sptr sceneA,boxm2_scene_sptr sceneB,const bocl_device_sptr & device,int nbins,bool do_vary_scale)11 boxm2_ocl_hierarchical_reg::boxm2_ocl_hierarchical_reg( boxm2_opencl_cache_sptr & cache,
12 boxm2_scene_sptr sceneA,
13 boxm2_scene_sptr sceneB,
14 const bocl_device_sptr& device,
15 int nbins, bool do_vary_scale): boxm2_ocl_reg_mutual_info(cache,sceneA,sceneB,device,nbins, 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_reg::init(vnl_vector<double> const& mu, vnl_vector<double> const & cov)
21 {
22 mu_ = mu ;
23 cov_= cov;
24 mu_cost_ = this->mutual_info(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_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 (unsigned int level =0; level <3; 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.flush();
50 int searchwidth = numsamples_search_width[level];
51 vnl_vector<double> best_sample = iter->second;
52 int numsamples_per_best_particle = (int) std::pow((float)searchwidth,(float)params_to_vary);
53 for(size_t sampleno = 0 ; sampleno < numsamples_per_best_particle; sampleno++)
54 {
55 int cont=sampleno;
56 vnl_vector<double> curr_sample = best_sample;
57 for(size_t k = 0; k <params_to_vary; k++)
58 {
59 int t = params_to_vary - k - 1;
60 unsigned int var = cont/(unsigned int)std::pow((float)searchwidth,(float)t); // quotient
61 cont = cont - (unsigned int)std::pow((float)searchwidth,(float)t)*var; // remainder
62 double offset = (2*((double)var/((double)searchwidth-1))-1)*cov[t];
63 curr_sample[t] += offset;
64 }
65 double minfo = this->mutual_info(curr_sample, level) ;
66 mis.push_back(minfo);
67 samples_.push_back(curr_sample);
68
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->mutual_info(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_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