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