1 // This is brl/bbas/bpgl/ihog/ihog_minfo_cost_func.cxx
2 #include <utility>
3 #include "ihog_minfo_cost_func.h"
4 //:
5 // \file
6 
7 #include <vil/algo/vil_gauss_filter.h>
8 #include "vbl/vbl_array_2d.h"
9 #include "vbl/vbl_array_1d.h"
10 
11 
12 //: Constructor
ihog_minfo_cost_func(const ihog_image<float> & image1,const ihog_image<float> & image2,ihog_world_roi roi,const ihog_transform_2d & init_xform,unsigned nbins)13 ihog_minfo_cost_func::ihog_minfo_cost_func( const ihog_image<float>& image1,
14                                             const ihog_image<float>& image2,
15                                             ihog_world_roi  roi,
16                                             const ihog_transform_2d& init_xform,
17                                             unsigned nbins)
18  : //vnl_least_squares_function(1,1),
19    vnl_cost_function(2),
20    from_image_(image1),
21    to_image_(image2),
22    roi_(std::move(roi)),
23    form_(init_xform.form()),
24    from_mask_(false),
25    to_mask_(false),
26    nbins_(nbins)
27 {
28   vnl_vector<double> params;
29   init_xform.params(params);
30   from_samples_ = roi_.sample(from_image_);
31   //use_gradient_ = false;
32   //int number_of_residuals = from_samples_.size();
33 #if 0
34   int number_of_residuals = 1;  // just the mutual info
35   vnl_least_squares_function::init(params.size(), number_of_residuals);
36 #endif
37 }
38 
39 
ihog_minfo_cost_func(const ihog_image<float> & image1,const ihog_image<float> & image2,const ihog_image<float> & mask,ihog_world_roi roi,const ihog_transform_2d & init_xform,bool image1_mask,unsigned nbins)40 ihog_minfo_cost_func::ihog_minfo_cost_func( const ihog_image<float>& image1,
41                                             const ihog_image<float>& image2,
42                                             const ihog_image<float>& mask,
43                                             ihog_world_roi  roi,
44                                             const ihog_transform_2d& init_xform, bool image1_mask, unsigned nbins)
45  : //vnl_least_squares_function(1,1),
46    vnl_cost_function(2),
47    from_image_(image1),
48    to_image_(image2),
49    roi_(std::move(roi)),
50    form_(init_xform.form()),
51    from_mask_(image1_mask),
52    to_mask_(!image1_mask),
53    nbins_(nbins)
54 {
55   if (from_mask_) {
56     from_mask_image_ = mask;
57   }
58   else {
59     to_mask_image_ = mask;
60   }
61   vnl_vector<double> params;
62   init_xform.params(params);
63   from_samples_ = roi_.sample(from_image_);
64 #if 0
65   int number_of_residuals = from_samples_.size();
66       number_of_residuals = 1;  // just the mutual info
67   use_gradient_ = false;
68   vnl_least_squares_function::init(params.size(), number_of_residuals);
69 #endif
70 }
71 
ihog_minfo_cost_func(const ihog_image<float> & image1,const ihog_image<float> & image2,const ihog_image<float> & mask1,const ihog_image<float> & mask2,ihog_world_roi roi,const ihog_transform_2d & init_xform,unsigned nbins)72 ihog_minfo_cost_func::ihog_minfo_cost_func(const ihog_image<float>& image1,
73                                            const ihog_image<float>& image2,
74                                            const ihog_image<float>& mask1,
75                                            const ihog_image<float>& mask2,
76                                            ihog_world_roi  roi,
77                                            const ihog_transform_2d& init_xform, unsigned nbins)
78  : //vnl_least_squares_function(1,1),
79    vnl_cost_function(2),
80    from_image_(image1),
81    to_image_(image2),
82    from_mask_image_(mask1),
83    to_mask_image_(mask2),
84    roi_(std::move(roi)),
85    form_(init_xform.form()),
86    from_mask_(true),
87    to_mask_(true),
88    nbins_(nbins)
89 {
90   vnl_vector<double> params;
91   init_xform.params(params);
92   from_samples_ = roi_.sample(from_image_);
93 #if 0
94   int number_of_residuals = from_samples_.size();
95       number_of_residuals = 1;  // just the mutual info
96   use_gradient_ = false;
97   vnl_least_squares_function::init(params.size(), number_of_residuals);
98 #endif
99 }
100 
101 
102 //: The main function.
103 //  Given the parameter vector x, compute the vector of residuals fx.
104 //  Fx has been sized appropriately before the call.  it should have dimension 1
105 //void
106 //ihog_minfo_cost_func::f(vnl_vector<double> const& x, vnl_vector<double>& fx)
f(vnl_vector<double> const & x)107 double ihog_minfo_cost_func::f(vnl_vector<double> const& x)
108 {
109   ihog_transform_2d new_xform;
110   new_xform.set(x, form_);
111   ihog_image<float> test_image(to_image_);
112   test_image.set_world2im(new_xform*to_image_.world2im());
113   vnl_vector<double> to_samples = roi_.sample(test_image);
114 
115   if (from_mask_ || to_mask_) {
116     vnl_vector<double> mask_samples;
117     if (from_mask_) {
118       mask_samples = roi_.sample(from_mask_image_);
119     }
120     if (to_mask_) {
121       ihog_image<float> mask_test_image(to_mask_image_);
122       mask_test_image.set_world2im(new_xform*to_image_.world2im());
123       vnl_vector<double> to_mask_samples = roi_.sample(mask_test_image);
124       if (from_mask_) {
125         // mask_samples already filled in - multiply by second mask to get AND
126         mask_samples = element_product<double>(mask_samples,to_mask_samples);
127       }
128       else {
129         // no previous mask - mask is to_mask
130         mask_samples = to_mask_samples;
131       }
132     }
133     //fx[0] = entropy_diff(mask_samples, to_samples);
134     return entropy_diff(mask_samples, from_samples_, to_samples, nbins_);
135   }
136   else
137   {
138     vnl_vector<double> mask_samples(from_samples_.size());
139     mask_samples.fill(1.0f);
140     //fx[0] = entropy_diff(mask_samples, to_samples);
141     return entropy_diff(mask_samples, from_samples_, to_samples, nbins_);
142   }
143 }
144 
145 
146 //: Returns the transformed second image
147 vil_image_view<float>
last_xformed_image()148 ihog_minfo_cost_func::last_xformed_image()
149 {
150   return roi_.resample(to_image_);
151 }
152 
153 
entropy_diff(vnl_vector<double> & mask_samples,vnl_vector<double> & from_samples,vnl_vector<double> & to_samples,int nbins)154 double ihog_minfo_cost_func::entropy_diff(vnl_vector<double>& mask_samples, vnl_vector<double>& from_samples, vnl_vector<double>& to_samples, int nbins)
155 {
156   double scl = 1.0/(256.0/nbins);
157   vbl_array_2d<double> h(nbins, nbins, 0.0);
158 
159   //compute the intensity histogram
160   double total_weight = 0.0;
161   for (unsigned i = 0; i<to_samples.size(); ++i)
162     if (mask_samples[i]>0.0) {
163       //match the gpu implementation, which does a floor operation
164       auto id = static_cast<unsigned>(std::floor(from_samples[i]*scl)),
165                is = static_cast<unsigned>(std::floor(to_samples[i]*scl));
166 
167       if (id+1>(unsigned)nbins || is+1>(unsigned)nbins)
168         continue;
169       h[id][is] += 1.0;
170       total_weight += 1.0;
171     }
172   // convert to probability
173   for (int r = 0; r<nbins; ++r)
174     for (int c = 0; c<nbins; ++c)
175       h[r][c] /= total_weight;
176 
177   auto nr = (unsigned)h.rows(), nc = (unsigned)h.cols();
178   //marginal distribution for mapped dest intensities
179   vbl_array_1d<double> pmr(nc,0.0);
180   for (unsigned r = 0; r<nr; ++r)
181     for (unsigned c = 0; c<nc; ++c)
182       pmr[c]+=h[r][c];
183   double jsum = 0.0, msum = 0.0;
184   for (unsigned c = 0; c<nc; ++c)
185   {
186     double pr = pmr[c];
187     if (pr>0)
188       msum += pr*std::log(pr);
189   }
190   for (unsigned r = 0; r<nr; ++r)
191     for (unsigned c = 0; c<nc; ++c) {
192         double prc = h[r][c];
193         if (prc>0)
194           jsum+= prc*std::log(prc);
195     }
196   double ent_dif = jsum - msum;
197   return -ent_dif/std::log(2.0);
198 }
199