1 #include <fstream>
2 #include <iostream>
3 #include <complex>
4 #include <limits>
5 #include "brip_vil_float_ops.h"
6 //:
7 // \file
8
9 #include <cassert>
10 #ifdef _MSC_VER
11 # include "vcl_msvc_warnings.h"
12 #endif
13 #include "vgl/vgl_box_2d.h"
14 #include "vgl/vgl_point_2d.h"
15 #include "vgl/vgl_vector_2d.h"
16 #include "vgl/vgl_polygon_scan_iterator.h"
17 #include "vul/vul_timer.h"
18 #include "vbl/vbl_array_1d.h"
19 #include "vbl/vbl_array_2d.h"
20 #include "vbl/vbl_bounding_box.h"
21 #include "vnl/vnl_numeric_traits.h"
22 #include "vnl/vnl_math.h"
23 #include "vnl/vnl_double_2x3.h"
24 #include <vnl/algo/vnl_fft_prime_factors.h>
25 #include <vnl/algo/vnl_svd.h>
26
27 #include "vil/vil_pixel_format.h"
28 #include "vil/vil_transpose.h"
29 #include "vil/vil_convert.h"
30 #include "vil/vil_save.h"
31 #include "vil/vil_new.h"
32 #include "vil/vil_math.h"
33 #include <vil/algo/vil_convolve_1d.h>
34
35 #include <vsol/vsol_box_2d.h>
36 #include <vsol/vsol_polygon_2d_sptr.h>
37 #include <vsol/vsol_polygon_2d.h>
38 #include <bsol/bsol_algs.h>
39 #include <bsta/bsta_histogram.h>
40 #include <bsta/bsta_joint_histogram.h>
41 #include <brip/brip_roi.h>
42
43 // === Local utility functions ===
44
45 //: compute normalized cross correlation from the intensity moment sums.
cross_corr(double area,double si1,double si2,double si1i1,double si2i2,double si1i2,float intensity_thresh)46 static float cross_corr(double area, double si1, double si2,
47 double si1i1, double si2i2, double si1i2,
48 float intensity_thresh)
49 {
50 if (!area)
51 return 0.f;
52 // the mean values
53 double u1 = si1/area, u2 = si2/area;
54 if (u1<intensity_thresh||u2<intensity_thresh)
55 return -1.f;
56 double neu = si1i2 - area*u1*u2;
57 double sd1 = std::sqrt(std::fabs(si1i1-area*u1*u1)),
58 sd2 = std::sqrt(std::fabs(si2i2-area*u2*u2));
59 if (!neu)
60 return 0.f;
61 if (!sd1||!sd2) {
62 if (neu>0)
63 return 1.f;
64 else
65 return -1.f;
66 }
67 double den = sd1*sd2;
68 return float(neu/den);
69 }
70
71 //------------------------------------------------------------
72 //: Convolve with a kernel
73 // It's assumed that the kernel is square with odd dimensions
74 vil_image_view<float>
convolve(vil_image_view<float> const & input,vbl_array_2d<float> const & kernel)75 brip_vil_float_ops::convolve(vil_image_view<float> const& input,
76 vbl_array_2d<float> const& kernel)
77 {
78 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
79 int kw = kernel.cols(); // kh = kernel.rows();
80 // add a check for kernels that are not equal dimensions of odd size JLM
81 int n = (kw-1)/2;
82 vil_image_view<float> output;
83 output.set_size(w,h);
84 for (int y = n; y<(h-n); y++)
85 for (int x = n; x<(w-n); x++)
86 {
87 float accum = 0;
88 for (int j = -n; j<=n; j++)
89 for (int i = -n; i<=n; i++)
90 {
91 float x1 = input(x+i,y+j);
92 float x2 = kernel[i+n][j+n];
93 accum += x1*x2;
94 }
95 output(x,y)=accum;
96 }
97 brip_vil_float_ops::fill_x_border(output, n, 0.0f);
98 brip_vil_float_ops::fill_y_border(output, n, 0.0f);
99 return output;
100 }
101
fill_1d_array(vil_image_view<float> const & input,int y,float * output)102 static void fill_1d_array(vil_image_view<float> const& input,
103 int y, float* output)
104 {
105 unsigned w = input.ni();
106 for (unsigned x = 0; x<w; x++)
107 output[x] = input(x,y);
108 }
109
110 //: Downsamples the 1-d array by 2 using the Burt-Adelson reduction algorithm.
half_resolution_1d(const float * input,unsigned width,float k0,float k1,float k2,float * output)111 void brip_vil_float_ops::half_resolution_1d(const float* input, unsigned width,
112 float k0, float k1, float k2,
113 float* output)
114 {
115 float w[5];
116 int n = 0;
117 for (; n<5; n++)
118 w[n]=input[n];
119 output[0]=k0*w[0]+ 2.0f*(k1*w[1] + k2*w[2]); // reflect at boundary
120 for (unsigned x = 1; x<width; ++x)
121 {
122 output[x]=k0*w[2]+ k1*(w[1]+w[3]) + k2*(w[0]+w[4]);
123 //shift the window, w, over by two pixels
124 w[0] = w[2]; w[1] = w[3]; w[2] = w[4];
125 //handle the boundary conditions
126 if (x+2<width)
127 w[3] = input[n++], w[4] = input[n++];
128 else
129 w[3] = w[1], w[4] = w[0];
130 }
131 }
132
133 //: Downsamples the image by 2 using the Burt-Adelson reduction algorithm.
134 // Convolution with a 5-point kernel [(0.5-ka)/2, 0.25, ka, 0.25, (0.5-ka)/2]
135 // ka = 0.6 maximum decorrelation, wavelet for image compression.
136 // ka = 0.5 linear interpolation,
137 // ka = 0.4 Gaussian filter
138 // ka = 0.359375 min aliasing, wider than Gaussian
139 // The image sizes are related by: output_dimension = (input_dimension +1)/2.
140 vil_image_view<float>
half_resolution(vil_image_view<float> const & input,float filter_coef)141 brip_vil_float_ops::half_resolution(vil_image_view<float> const& input,
142 float filter_coef)
143 {
144 vul_timer t;
145 float k0 = filter_coef, k1 = 0.25f*filter_coef, k2 = 0.5f*(0.5f-filter_coef);
146 unsigned w = input.ni(), h = input.nj();
147 int half_w =(w+1)/2, half_h = (h+1)/2;
148 vil_image_view<float> output;
149 output.set_size(half_w, half_h);
150 // Generate input/output arrays
151 int n = 0;
152 auto* in0 = new float[w]; auto* in1 = new float[w];
153 auto* in2 = new float[w]; auto* in3 = new float[w];
154 auto* in4 = new float[w];
155
156 auto* out0 = new float[half_w]; auto* out1 = new float[half_w];
157 auto* out2 = new float[half_w]; auto* out3 = new float[half_w];
158 auto* out4 = new float[half_w];
159 // Initialize arrays
160 fill_1d_array(input, n++, in0); fill_1d_array(input, n++, in1);
161 fill_1d_array(input, n++, in2); fill_1d_array(input, n++, in3);
162 fill_1d_array(input, n++, in4);
163
164 // downsample initial arrays
165 brip_vil_float_ops::half_resolution_1d(in0, half_w, k0, k1, k2, out0);
166 brip_vil_float_ops::half_resolution_1d(in1, half_w, k0, k1, k2, out1);
167 brip_vil_float_ops::half_resolution_1d(in2, half_w, k0, k1, k2, out2);
168 brip_vil_float_ops::half_resolution_1d(in3, half_w, k0, k1, k2, out3);
169 brip_vil_float_ops::half_resolution_1d(in4, half_w, k0, k1, k2, out4);
170 int x=0, y;
171 // do the first output line
172 for (;x<half_w;x++)
173 output(x,0)= k0*out0[x]+ 2.0f*(k1*out1[x]+k2*out2[x]);
174 // normal lines
175 for (y=1; y<half_h; y++)
176 {
177 for (x=0; x<half_w; x++)
178 output(x,y) = k0*out2[x]+ k1*(out1[x]+out3[x]) + k2*(out0[x]+out4[x]);
179 //shift the neighborhood down two lines
180 float* temp0 = out0;
181 float* temp1 = out1;
182 out0 = out2; out1 = out3; out2 = out4;
183 out3 = temp0; out4 = temp1;//reflect values
184 //test border condition
185 if (y<half_h-2)
186 {
187 //normal processing, so don't reflect
188 fill_1d_array(input, n++, in3);
189 fill_1d_array(input, n++, in4);
190 brip_vil_float_ops::half_resolution_1d(in3, half_w, k0, k1, k2, out3);
191 brip_vil_float_ops::half_resolution_1d(in4, half_w, k0, k1, k2, out4);
192 }
193 }
194 delete [] in0; delete [] in1; delete [] in2;
195 delete [] in3; delete [] in4;
196 delete [] out0; delete [] out1; delete [] out2;
197 delete [] out3; delete [] out4;
198 #ifdef DEBUG
199 std::cout << "\nDownsample a "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
200 #endif
201 return output;
202 }
203
double_resolution_1d(const float * input,const unsigned n_input,const float k0,const float k1,const float k2,float * output)204 void brip_vil_float_ops::double_resolution_1d(const float* input, const unsigned n_input,
205 const float k0, const float k1,
206 const float k2, float* output)
207 {
208 float w[3];
209 unsigned i = 0;
210 w[1]=input[i]; w[2]=input[i++];
211 w[0]=w[2];
212 for (unsigned c = 0; c<2*n_input; c+=2)
213 {
214 output[c] = k0*w[1] + k2*(w[0]+w[2]);
215 output[c+1] = k1*(w[1]+w[2]);
216 w[0]=w[1];
217 w[1]=w[2];
218 if (c<2*(n_input-2))
219 w[2]=input[i++];
220 else
221 w[2]=w[0];
222 }
223 }
224
225 //: interpolates the input using the Bert-Adelson algorithm
226 vil_image_view<float>
double_resolution(vil_image_view<float> const & input,float filter_coef)227 brip_vil_float_ops::double_resolution(vil_image_view<float> const& input,
228 float filter_coef)
229 {
230 unsigned ni_in = input.ni();
231 unsigned nj_in = input.nj();
232 unsigned ni_out = 2*ni_in;
233 unsigned nj_out = 2*nj_in;
234 vil_image_view<float> out(ni_out, nj_out);
235 auto* input_1d = new float[ni_in];
236
237 // An interpolation neighborhood of three lines
238 auto* output0 = new float[ni_out];
239 auto* output1 = new float[ni_out];
240 auto* output2 = new float[ni_out];
241
242 // The filter coefficients
243 float k0 = filter_coef*2.0f;
244 float k1 = 0.5f;
245 float k2 = 0.5f-filter_coef;
246
247 // initialize
248 unsigned i = 0;
249 fill_1d_array(input, i++, input_1d);
250 brip_vil_float_ops::double_resolution_1d(input_1d, ni_in, k0, k1, k2, output1);
251 fill_1d_array(input, i++, input_1d);
252 brip_vil_float_ops::double_resolution_1d(input_1d, ni_in, k0, k1, k2, output2);
253 for (unsigned k = 0; k<ni_out; ++k)
254 output0[k]=output2[k];
255 for (unsigned r = 0; r<nj_out; r+=2)
256 {
257 unsigned rp = r+1;
258 for (unsigned c=0; c<ni_out; ++c)
259 {
260 out(c, r) = k0*output1[c] + k2*(output0[c]+output2[c]);
261 out(c, rp) = k1*(output1[c]+output2[c]);
262 }
263 float* next = output0;
264 output0 = output1;
265 output1 = output2;
266 output2 = next;
267 if (r<nj_out-4)
268 {
269 fill_1d_array(input, i++, input_1d);
270 brip_vil_float_ops::double_resolution_1d(input_1d, ni_in,
271 k0, k1, k2, output2);
272 }
273 else
274 for (unsigned k = 0; k<ni_out; ++k)
275 output2[k]=output0[k];
276 }
277 delete [] input_1d;
278 delete [] output0;
279 delete [] output1;
280 delete [] output2;
281 return out;
282 }
283
brip_vil_gaussian(double x,double sigma)284 static double brip_vil_gaussian(double x, double sigma)
285 {
286 double x_on_sigma = x / sigma;
287 return (double)std::exp(- x_on_sigma * x_on_sigma / 2);
288 }
289
gaussian_radius(const double sigma,const double fuzz)290 unsigned brip_vil_float_ops::gaussian_radius(const double sigma,
291 const double fuzz)
292 {
293 unsigned radius = 0;
294 while ((float)brip_vil_gaussian((double)radius, sigma) > fuzz) ++radius;
295 return radius;
296 }
297
298 //: generate a 1-d Gaussian kernel fuzz=0.02 is a good value
brip_1d_gaussian_kernel(double sigma,double fuzz,unsigned & radius,double * & kernel)299 static void brip_1d_gaussian_kernel(double sigma, double fuzz,
300 unsigned& radius, double*& kernel)
301 {
302 radius = brip_vil_float_ops::gaussian_radius(sigma, fuzz);
303
304 kernel = new double[2*radius + 1];
305 if (!radius)
306 {
307 kernel[0]=1;
308 return;
309 }
310 for (unsigned i=0; i<=radius; ++i)
311 kernel[radius+i] = kernel[radius-i] = brip_vil_gaussian(double(i), sigma);
312 double sum = 0;
313 for (unsigned i= 0; i <= 2*radius; ++i)
314 sum += kernel[i]; // find integral of weights
315 for (unsigned i= 0; i <= 2*radius; ++i)
316 kernel[i] /= sum; // normalize by integral
317 }
318
319 vil_image_view<float>
gaussian(vil_image_view<float> const & input,float sigma,std::string const & boundary_condition,float fill)320 brip_vil_float_ops::gaussian(vil_image_view<float> const& input, float sigma,
321 std::string const& boundary_condition,
322 float fill)
323 {
324 vil_convolve_boundary_option option=vil_convolve_ignore_edge;
325 if(boundary_condition == "zeros")
326 option = vil_convolve_zero_extend;
327 else if(boundary_condition == "const")
328 option = vil_convolve_constant_extend;
329 else if(boundary_condition == "periodic")
330 option = vil_convolve_periodic_extend;
331 else if(boundary_condition == "reflect")
332 option = vil_convolve_reflect_extend;
333
334 unsigned ni = input.ni(), nj = input.nj();
335 unsigned np = input.nplanes();
336 vil_image_view<float> dest(ni, nj, np);
337 dest.fill(fill);
338 unsigned r;
339 double* ker;
340 brip_1d_gaussian_kernel(sigma, 0.02, r, ker);
341 for (unsigned p = 0; p<np; ++p) {
342 vil_image_view<float> input_temp(ni, nj);
343 vil_image_view<float> out_temp(ni, nj);
344 vil_image_view<float> work(ni, nj);
345 work.fill(fill);
346 for (unsigned j = 0; j<nj; ++j)
347 for (unsigned i = 0; i<ni; ++i)
348 input_temp(i,j) = input(i,j,p);
349
350 // filter horizontal
351 int ksize = 2*r + 1 ;
352 float accum=0.0f;
353 vil_convolve_1d(input_temp, work, ker + ksize/2,
354 -ksize/2, r, accum,
355 option, option);
356 // filter vertical
357 vil_image_view<float> work_t = vil_transpose(work);
358 vil_image_view<float> out_temp_t = vil_transpose(dest);
359 vil_convolve_1d(work_t, out_temp_t, ker+ ksize/2,
360 -ksize/2, r, accum,
361 option, option);
362
363 vil_image_view<float> plane = vil_transpose(out_temp_t);
364 for (unsigned j = 0; j<nj; ++j)
365 for (unsigned i = 0; i<ni; ++i)
366 dest(i,j,p) = plane(i,j);
367 }
368 delete [] ker;
369 return dest;
370 }
371
372 vil_image_view<float>
absolute_value(vil_image_view<float> const & input)373 brip_vil_float_ops::absolute_value(vil_image_view<float> const& input)
374 {
375 unsigned ni = input.ni(), nj = input.nj();
376 unsigned np = input.nplanes();
377 vil_image_view<float> dest(ni, nj, np);
378 for (unsigned p = 0; p<np; ++p)
379 for (unsigned j = 0; j<nj; ++j)
380 for (unsigned i = 0; i<ni; ++i)
381 dest(i,j,p) = std::fabs(input(i,j,p));
382 return dest;
383 }
384
385 #ifdef VIL_CONVOLVE_WITH_MASK_EXISTS // TODO
386 vil_image_view<float>
gaussian(vil_image_view<float> const & input,float sigma,vil_image_view<float> const & mask)387 brip_vil_float_ops::gaussian(vil_image_view<float> const& input,
388 float sigma,
389 vil_image_view<float> const& mask)
390 {
391 vil_image_view<float> dest(input.ni(), input.nj());
392
393 int r;
394 double* ker;
395 brip_1d_gaussian_kernel(sigma, 0.02, r, ker);
396 vil_image_view<float> work(input.ni(), input.nj());
397 work.deep_copy(input);
398 // filter horizontal
399 int ksize = 2*r + 1 ;
400 float accum=0.0f;
401 vil_convolve_1d(input, work, mask, ker + ksize/2,
402 -ksize/2, r, accum,
403 vil_convolve_trim,
404 vil_convolve_trim);
405
406 // filter vertical
407 vil_image_view<float> work_t = vil_transpose(work);
408 vil_image_view<float> dest_t = vil_transpose(dest);
409 vil_convolve_1d(work_t, dest_t, vil_transpose(mask), ker+ ksize/2,
410 -ksize/2, r, accum,
411 vil_convolve_constant_extend,
412 vil_convolve_constant_extend);
413
414 delete ker;
415 return dest;
416 }
417
418 #endif // VIL_CONVOLVE_WITH_MASK_EXISTS
419
420 //-------------------------------------------------------------------
421 //: Determine if the center of a (2n+1)x(2n+1) neighborhood is a local maximum
422 //
423 bool brip_vil_float_ops::
local_maximum(vbl_array_2d<float> const & neighborhood,int n,float & value)424 local_maximum(vbl_array_2d<float> const& neighborhood, int n, float& value)
425 {
426 bool local_max = true;
427 value = 0;
428 float center = neighborhood[n][n];
429 for (int y = -n; y<=n; y++)
430 for (int x = -n; x<=n; x++)
431 local_max = local_max&&(neighborhood[y+n][x+n]<=center);
432 if (!local_max)
433 return false;
434 value = center;
435 return true;
436 }
437
438 //-------------------------------------------------------------------
439 // Interpolate the sub-pixel position of a neighborhood using a
440 // second order expansion on a 3x3 sub-neighborhood. Return the
441 // offset to the maximum, i.e. x=x0+dx, y = y0+dy. The design is
442 // similar to the droid counterpoint by fsm, which uses the Beaudet Hessian
443 //
444 void brip_vil_float_ops::
interpolate_center(vbl_array_2d<float> const & neighborhood,float & dx,float & dy)445 interpolate_center(vbl_array_2d<float>const& neighborhood, float& dx, float& dy)
446 {
447 dx = 0; dy=0;
448 // extract the neighborhood
449 float n_m1_m1 = neighborhood[0][0];
450 float n_m1_0 = neighborhood[0][1];
451 float n_m1_1 = neighborhood[0][2];
452 float n_0_m1 = neighborhood[1][0];
453 float n_0_0 = neighborhood[1][1];
454 float n_0_1 = neighborhood[1][2];
455 float n_1_m1 = neighborhood[2][0];
456 float n_1_0 = neighborhood[2][1];
457 float n_1_1 = neighborhood[2][2];
458
459 // Compute the 2nd order quadratic coefficients
460 // 1/6 * [ -1 0 +1 ]
461 // Ix = [ -1 0 +1 ]
462 // [ -1 0 +1 ]
463 float Ix =(-n_m1_m1+n_m1_1-n_0_m1+n_0_1-n_1_m1+n_1_1)/6.0f;
464 // 1/6 * [ -1 -1 -1 ]
465 // Iy = [ 0 0 0 ]
466 // [ +1 +1 +1 ]
467 float Iy =(-n_m1_m1-n_m1_0-n_m1_1+n_1_m1+n_1_0+n_1_1)/6.0f;
468 // 1/3 * [ +1 -2 +1 ]
469 // Ixx = [ +1 -2 +1 ]
470 // [ +1 -2 +1 ]
471 float Ixx = ((n_m1_m1+n_0_m1+n_1_m1+n_m1_1+n_0_1+n_1_1)
472 -2.0f*(n_m1_0+n_0_0+n_1_0))/3.0f;
473 // 1/4 * [ +1 0 -1 ]
474 // Ixy = [ 0 0 0 ]
475 // [ -1 0 +1 ]
476 float Ixy = (n_m1_m1-n_m1_1+n_1_m1+n_1_1)/4.0f;
477 // 1/3 * [ +1 +1 +1 ]
478 // Iyy = [ -2 -2 -2 ]
479 // [ +1 +1 +1 ]
480 float Iyy = ((n_m1_m1+n_m1_0+n_m1_1+n_1_m1+n_1_0+n_1_1)
481 -2.0f*(n_0_m1 + n_0_0 + n_1_0))/3.0f;
482 //
483 // The next bit is to find the extremum of the fitted surface by setting its
484 // partial derivatives to zero. We need to solve the following linear system :
485 // Given the fitted surface is
486 // I(x,y) = Io + Ix x + Iy y + 1/2 Ixx x^2 + Ixy x y + 1/2 Iyy y^2
487 // we solve for the maximum (x,y),
488 //
489 // [ Ixx Ixy ] [ dx ] + [ Ix ] = [ 0 ] (dI/dx = 0)
490 // [ Ixy Iyy ] [ dy ] [ Iy ] [ 0 ] (dI/dy = 0)
491 //
492 float det = Ixx*Iyy - Ixy*Ixy;
493 // det>0 corresponds to a true local extremum otherwise a saddle point
494 if (det>0)
495 {
496 dx = (Iy*Ixy - Ix*Iyy) / det;
497 dy = (Ix*Ixy - Iy*Ixx) / det;
498 // more than one pixel away
499 if (std::fabs(dx) > 1.0 || std::fabs(dy) > 1.0)
500 dx = 0; dy = 0;
501 }
502 }
503
504 //---------------------------------------------------------------
505 // Compute the local maxima of the input on a (2n+1)x(2n+2)
506 // neighborhood above the given threshold. At each local maximum,
507 // compute the sub-pixel location, (x_pos, y_pos).
508 void brip_vil_float_ops::
non_maximum_suppression(vil_image_view<float> const & input,int n,float thresh,std::vector<float> & x_pos,std::vector<float> & y_pos,std::vector<float> & value)509 non_maximum_suppression(vil_image_view<float> const& input,
510 int n, float thresh,
511 std::vector<float>& x_pos,
512 std::vector<float>& y_pos,
513 std::vector<float>& value)
514 {
515 vul_timer t;
516 int N = 2*n+1;
517 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
518 x_pos.clear(); x_pos.clear(); value.clear();
519 vbl_array_2d<float> neighborhood(N,N);
520 for (int y =n; y<h-n; y++)
521 for (int x = n; x<w-n; x++)
522 {
523 //If the center is not above threshold then there is
524 //no hope
525 if (input(x,y)<thresh)
526 continue;
527 //Fill the neighborhood
528 for (int i = -n; i<=n; i++)
529 for (int j = -n; j<=n; j++)
530 neighborhood.put(j+n,i+n,input(x+i, y+j));
531 //Check if the center is a local maximum
532 float dx, dy, max_v;
533 if (brip_vil_float_ops::local_maximum(neighborhood, n, max_v))
534 {
535 //if so sub-pixel interpolate (3x3) and output results
536 brip_vil_float_ops::interpolate_center(neighborhood, dx, dy);
537 x_pos.push_back((float)x+dx);
538 y_pos.push_back((float)y+dy);
539 value.push_back(max_v);
540 }
541 }
542 #ifdef DEBUG
543 std::cout << "\nCompute non-maximum suppression on a "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
544 #endif
545 }
546
547 // -----------------------------------------------------------------
548 //: Subtract image_1 from image_2.
549 // Will not operate unless the two input images are the same dimensions
550 //
551 vil_image_view<float>
difference(vil_image_view<float> const & image_1,vil_image_view<float> const & image_2)552 brip_vil_float_ops::difference(vil_image_view<float> const& image_1,
553 vil_image_view<float> const& image_2)
554 {
555 unsigned w1 = image_1.ni(), h1 = image_1.nj();
556 unsigned w2 = image_2.ni(), h2 = image_2.nj();
557 vil_image_view<float> temp(w1, h1);
558 if (w1!=w2||h1!=h2)
559 {
560 std::cout << "In brip_vil_float_ops::difference(..) - images are not the same dimensions\n";
561 return temp;
562 }
563 vil_image_view<float> out;
564 out.set_size(w1, h1);
565 for (unsigned y = 0; y<h1; y++)
566 for (unsigned x = 0; x<w1; x++)
567 out(x,y) = image_2(x,y)-image_1(x,y);
568 return out;
569 }
570
571 //: negate an image returning the same pixel type (only greyscale)
negate(vil_image_resource_sptr const & imgr)572 vil_image_resource_sptr brip_vil_float_ops::negate(vil_image_resource_sptr const& imgr)
573 {
574 vil_image_resource_sptr outr;
575 if (!imgr)
576 return outr;
577
578 vil_pixel_format fmt = imgr->pixel_format();
579 switch (fmt)
580 {
581 #define NEGATE_CASE(FORMAT, T) \
582 case FORMAT: { \
583 vil_image_view<T> view = imgr->get_copy_view(); \
584 T mxv = std::numeric_limits<T>::max(); \
585 vil_math_scale_and_offset_values(view, -1.0, mxv); \
586 outr = vil_new_image_resource_of_view(view); \
587 break; \
588 }
589 NEGATE_CASE(VIL_PIXEL_FORMAT_BYTE, vxl_byte);
590 NEGATE_CASE(VIL_PIXEL_FORMAT_UINT_32, vxl_uint_32);
591 NEGATE_CASE(VIL_PIXEL_FORMAT_UINT_16, vxl_uint_16);
592 NEGATE_CASE(VIL_PIXEL_FORMAT_INT_16, vxl_int_16);
593 NEGATE_CASE(VIL_PIXEL_FORMAT_FLOAT, float);
594 NEGATE_CASE(VIL_PIXEL_FORMAT_DOUBLE, double);
595 #undef NEGATE_CASE
596 default:
597 std::cout << "Unknown image format\n";
598 }
599 return outr;
600 }
601
602 vil_image_view<float>
threshold(vil_image_view<float> const & image,const float thresh,const float level)603 brip_vil_float_ops::threshold(vil_image_view<float> const & image,
604 const float thresh, const float level)
605 {
606 vil_image_view<float> out;
607 unsigned w = image.ni(), h = image.nj();
608 out.set_size(w, h);
609 for (unsigned y = 0; y<h; y++)
610 for (unsigned x = 0; x<w; x++)
611 {
612 if (image(x,y)>thresh)
613 out(x,y) = level;
614 else
615 out(x,y) = 0;
616 }
617 return out;
618 }
619
620 vil_image_view<float>
abs_clip_to_level(vil_image_view<float> const & image,float thresh,float level)621 brip_vil_float_ops::abs_clip_to_level(vil_image_view<float> const& image,
622 float thresh, float level)
623 {
624 vil_image_view<float> out;
625 unsigned w = image.ni(), h = image.nj();
626 out.set_size(w, h);
627 for (unsigned y = 0; y<h; y++)
628 for (unsigned x = 0; x<w; x++)
629 {
630 if (std::fabs(image(x,y))>thresh)
631 out(x,y) = level;
632 else
633 out(x,y) = image(x,y);
634 }
635 return out;
636 }
637
average_NxN(vil_image_view<float> const & img,int N)638 vil_image_view<float> brip_vil_float_ops::average_NxN(vil_image_view<float> const & img, int N)
639 {
640 vil_image_view<float> result;
641 unsigned w = img.ni(), h = img.nj();
642 result.set_size (w, h);
643
644 vbl_array_2d <float> averaging_filt(N, N);
645 averaging_filt.fill( float(1.00/double(N*N)) );
646 result = brip_vil_float_ops::convolve(img, averaging_filt);
647 return result;
648 }
649
650 //----------------------------------------------------------------
651 //: Compute the gradient of the input, use a 3x3 mask
652 // \verbatim
653 // 1 |-1 0 1| 1 |-1 -1 -1|
654 // Ix = --- |-1 0 1| Iy = --- | 0 0 0|
655 // 6 |-1 0 1| 6 | 1 1 1|
656 // \endverbatim
657 // Larger masks are computed by pre-convolving with a Gaussian
658 //
gradient_3x3(vil_image_view<float> const & input,vil_image_view<float> & grad_x,vil_image_view<float> & grad_y)659 void brip_vil_float_ops::gradient_3x3(vil_image_view<float> const& input,
660 vil_image_view<float>& grad_x,
661 vil_image_view<float>& grad_y)
662 {
663 vul_timer t;
664 unsigned w = input.ni(), h = input.nj();
665 float scale = 1.0f/6.0f;
666 for (unsigned y = 1; y+1<h; ++y)
667 for (unsigned x = 1; x+1<w; ++x)
668 {
669 float gx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
670 -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
671 float gy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
672 -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
673 grad_x(x,y) = scale*gx;
674 grad_y(x,y) = scale*gy;
675 }
676 brip_vil_float_ops::fill_x_border(grad_x, 1, 0.0f);
677 brip_vil_float_ops::fill_y_border(grad_x, 1, 0.0f);
678 brip_vil_float_ops::fill_x_border(grad_y, 1, 0.0f);
679 brip_vil_float_ops::fill_y_border(grad_y, 1, 0.0f);
680 #ifdef DEBUG
681 std::cout << "\nCompute Gradient in " << t.real() << " msecs.\n";
682 #endif
683 }
684
gradient_mag_3x3(vil_image_view<float> const & input,vil_image_view<float> & mag)685 void brip_vil_float_ops::gradient_mag_3x3(vil_image_view<float> const& input,
686 vil_image_view<float>& mag)
687 {
688 unsigned w = input.ni(), h = input.nj();
689 float scale = 1.0f/6.0f;
690 for (unsigned y = 1; y+1<h; ++y)
691 for (unsigned x = 1; x+1<w; ++x)
692 {
693 float gx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
694 -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
695 float gy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
696 -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
697 mag(x,y) = scale*std::sqrt(gx*gx+gy*gy);
698 }
699 brip_vil_float_ops::fill_x_border(mag, 1, 0.0f);
700 brip_vil_float_ops::fill_y_border(mag, 1, 0.0f);
701 }
702
703 void brip_vil_float_ops::
gradient_mag_comp_3x3(vil_image_view<float> const & input,vil_image_view<float> & mag,vil_image_view<float> & gx,vil_image_view<float> & gy)704 gradient_mag_comp_3x3(vil_image_view<float> const& input,
705 vil_image_view<float>& mag,
706 vil_image_view<float>& gx,
707 vil_image_view<float>& gy)
708 {
709 unsigned w = input.ni(), h = input.nj();
710 float scale = 1.0f/6.0f;
711 for (unsigned y = 1; y+1<h; ++y)
712 for (unsigned x = 1; x+1<w; ++x)
713 {
714 float ggx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
715 -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
716 float ggy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
717 -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
718 mag(x,y) = scale*std::sqrt(ggx*ggx+ggy*ggy);
719 gx(x,y) = ggx*scale;
720 gy(x,y) = ggy*scale;
721 }
722 brip_vil_float_ops::fill_x_border(mag, 1, 0.0f);
723 brip_vil_float_ops::fill_y_border(mag, 1, 0.0f);
724 brip_vil_float_ops::fill_x_border(gx, 1, 0.0f);
725 brip_vil_float_ops::fill_y_border(gy, 1, 0.0f);
726 }
727
728 //----------------------------------------------------------------
729 //: Compute the Hessian of the input, use a 3x3 mask
730 // \verbatim
731 // 1 | 1 -2 1| 1 | 1 1 1| 1 | 1 0 -1|
732 // Ixx = ---| 1 -2 1| Iyy = ---| -2 -2 -2| Ixy = --- | 0 0 0|
733 // 3 | 1 -2 1| 3 | 1 1 1| 4 |-1 0 1|
734 // \endverbatim
735 // Larger masks are computed by pre-convolving with a Gaussian
736 //
hessian_3x3(vil_image_view<float> const & input,vil_image_view<float> & Ixx,vil_image_view<float> & Ixy,vil_image_view<float> & Iyy)737 void brip_vil_float_ops::hessian_3x3(vil_image_view<float> const& input,
738 vil_image_view<float>& Ixx,
739 vil_image_view<float>& Ixy,
740 vil_image_view<float>& Iyy)
741 {
742 vul_timer t;
743 unsigned w = input.ni(), h = input.nj();
744 for (unsigned y = 1; y+1<h; ++y)
745 for (unsigned x = 1; x+1<w; ++x)
746 {
747 float xx = input(x-1,y-1)+input(x-1,y)+input(x+1,y+1)+
748 input(x+1,y-1)+input(x+1,y)+input(x+1,y+1)-
749 2.0f*(input(x,y-1)+input(x,y)+input(x,y+1));
750
751 float xy = (input(x-1,y-1)+input(x+1,y+1))-
752 (input(x-1,y+1)+input(x+1,y-1));
753
754 float yy = input(x-1,y-1)+input(x,y-1)+input(x+1,y-1)+
755 input(x-1,y+1)+input(x,y+1)+input(x+1,y+1)-
756 2.0f*(input(x-1,y)+input(x,y)+input(x+1,y));
757
758 Ixx(x,y) = xx/3.0f;
759 Ixy(x,y) = xy/4.0f;
760 Iyy(x,y) = yy/3.0f;
761 }
762 brip_vil_float_ops::fill_x_border(Ixx, 1, 0.0f);
763 brip_vil_float_ops::fill_y_border(Ixx, 1, 0.0f);
764 brip_vil_float_ops::fill_x_border(Ixy, 1, 0.0f);
765 brip_vil_float_ops::fill_y_border(Ixy, 1, 0.0f);
766 brip_vil_float_ops::fill_x_border(Iyy, 1, 0.0f);
767 brip_vil_float_ops::fill_y_border(Iyy, 1, 0.0f);
768 #ifdef DEBUG
769 std::cout << "\nCompute a hessian matrix "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
770 #endif
771 }
772
773 vil_image_view<float>
beaudet(vil_image_view<float> const & Ixx,vil_image_view<float> const & Ixy,vil_image_view<float> const & Iyy,bool determinant)774 brip_vil_float_ops::beaudet(vil_image_view<float> const& Ixx,
775 vil_image_view<float> const& Ixy,
776 vil_image_view<float> const& Iyy,
777 bool determinant)
778 {
779 unsigned w = Ixx.ni(), h = Ixx.nj();
780 vil_image_view<float> output;
781 output.set_size(w, h);
782 for (unsigned y = 0; y<h; y++)
783 for (unsigned x = 0; x<w; x++)
784 {
785 float xx = Ixx(x,y), xy = Ixy(x,y), yy = Iyy(x,y);
786
787 //compute eigenvalues for experimentation
788 float det = xx*yy-xy*xy;
789 float tr = xx+yy;
790 if (determinant)
791 output(x,y) = det;
792 else
793 output(x,y) = tr;
794 }
795 return output;
796 }
797
798 //----------------------------------------------------------------
799 //: $Ix\cdot Ix^t$ gradient matrix elements
800 // That is,
801 // \verbatim
802 // _ _
803 // | (dI/dx)^2 (dI/dx)(dI/dy) |
804 // | |
805 // A = Sum(neighborhood)| |
806 // |(dI/dx)(dI/dy) (dI/dx)^2 |
807 // |_ _|
808 // \endverbatim
809 // over a 2n+1 x 2n+1 neighborhood
810 //
811 void
grad_matrix_NxN(vil_image_view<float> const & input,unsigned n,vil_image_view<float> & IxIx,vil_image_view<float> & IxIy,vil_image_view<float> & IyIy)812 brip_vil_float_ops::grad_matrix_NxN(vil_image_view<float> const& input,
813 unsigned n,
814 vil_image_view<float>& IxIx,
815 vil_image_view<float>& IxIy,
816 vil_image_view<float>& IyIy)
817 {
818 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
819 int N = (2*n+1)*(2*n+1);
820 int ni = static_cast<int>(n);
821 vil_image_view<float> grad_x, grad_y, output;
822 grad_x.set_size(w,h);
823 grad_y.set_size(w,h);
824 output.set_size(w,h);
825 brip_vil_float_ops::gradient_3x3(input, grad_x, grad_y);
826 vul_timer t;
827 for (int y = ni; y<h-ni;y++)
828 for (int x = ni; x<w-ni;x++)
829 {
830 float xx=0, xy=0, yy=0;
831 for (int i = -ni; i<=ni; i++)
832 for (int j = -ni; j<=ni; j++)
833 {
834 float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
835 xx += gx*gx;
836 xy += gx*gy;
837 yy += gy*gy;
838 }
839 IxIx(x,y) = xx/(float)N;
840 IxIy(x,y) = xy/(float)N;
841 IyIy(x,y) = yy/(float)N;
842 }
843 brip_vil_float_ops::fill_x_border(IxIx, ni, 0.0f);
844 brip_vil_float_ops::fill_y_border(IxIx, ni, 0.0f);
845 brip_vil_float_ops::fill_x_border(IxIy, ni, 0.0f);
846 brip_vil_float_ops::fill_y_border(IxIy, ni, 0.0f);
847 brip_vil_float_ops::fill_x_border(IyIy, ni, 0.0f);
848 brip_vil_float_ops::fill_y_border(IyIy, ni, 0.0f);
849 #ifdef DEBUG
850 std::cout << "\nCompute a gradient matrix "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
851 #endif
852 }
853
854 vil_image_view<float> brip_vil_float_ops::
trace_grad_matrix_NxN(vil_image_view<float> const & input,unsigned n)855 trace_grad_matrix_NxN(vil_image_view<float> const& input, unsigned n)
856 {
857 unsigned ni = input.ni(), nj = input.nj();
858 vil_image_view<float> IxIx;
859 vil_image_view<float> IxIy;
860 vil_image_view<float> IyIy;
861 vil_image_view<float> tr;
862 IxIx.set_size(ni, nj); IxIy.set_size(ni, nj); IyIy.set_size(ni, nj);
863 tr.set_size(ni, nj);
864 brip_vil_float_ops::grad_matrix_NxN(input, n, IxIx, IxIy, IyIy);
865 vil_math_image_sum<float, float, float>(IxIx, IyIy, tr);
866 return tr;
867 }
868
869 vil_image_view<float>
harris(vil_image_view<float> const & IxIx,vil_image_view<float> const & IxIy,vil_image_view<float> const & IyIy,double scale)870 brip_vil_float_ops::harris(vil_image_view<float> const& IxIx,
871 vil_image_view<float> const& IxIy,
872 vil_image_view<float> const& IyIy,
873 double scale)
874 {
875 unsigned w = IxIx.ni(), h = IxIx.nj();
876 float norm = 1e-3f; // Scale the output to values in the 10->1000 range
877 vil_image_view<float> output;
878 output.set_size(w, h);
879 for (unsigned y = 0; y<h; y++)
880 for (unsigned x = 0; x<w; x++)
881 {
882 float xx = IxIx(x,y), xy = IxIy(x,y), yy = IyIy(x,y);
883 float det = xx*yy-xy*xy, trace = xx+yy;
884 output(x,y) = float(det - scale*trace*trace)*norm;
885 }
886 return output;
887 }
888
889 //----------------------------------------------------------------
890 // Compute the sqrt of the product of the eigenvalues of the
891 // gradient matrix over a 2n+1 x 2n+1 neighborhood
892 // That is,
893 // \verbatim
894 // _ _
895 // | (dI/dx)^2 (dI/dx)(dI/dy) |
896 // | |
897 // A = Sum(neighborhood)| |
898 // |(dI/dx)(dI/dy) (dI/dx)^2 |
899 // |_ _|
900 // \endverbatim
901 // The output image is sqrt(lamba_1*lambda_2) where lambda_i are the eigenvalues
902 //
903 vil_image_view<float>
sqrt_grad_singular_values(vil_image_view<float> & input,int n)904 brip_vil_float_ops::sqrt_grad_singular_values(vil_image_view<float>& input,
905 int n)
906 {
907 int N = (2*n+1)*(2*n+1);
908 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
909 vil_image_view<float> grad_x, grad_y, output;
910 grad_x.set_size(w,h);
911 grad_y.set_size(w,h);
912 output.set_size(w,h);
913 brip_vil_float_ops::gradient_3x3(input, grad_x, grad_y);
914 vul_timer t;
915 for (int y = n; y<h-n;y++)
916 for (int x = n; x<w-n;x++)
917 {
918 float IxIx=0, IxIy=0, IyIy=0;
919 for (int i = -n; i<=n; i++)
920 for (int j = -n; j<=n; j++)
921 {
922 float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
923 IxIx += gx*gx;
924 IxIy += gx*gy;
925 IyIy += gy*gy;
926 }
927 float det = (IxIx*IyIy-IxIy*IxIy)/(float)N;
928 output(x,y)=std::sqrt(std::fabs(det));
929 }
930 brip_vil_float_ops::fill_x_border(output, n, 0.0f);
931 brip_vil_float_ops::fill_y_border(output, n, 0.0f);
932 #ifdef DEBUG
933 std::cout << "\nCompute sqrt(sigma0*sigma1) in" << t.real() << " msecs.\n";
934 #endif
935 return output;
936 }
937
938 vil_image_view<float> brip_vil_float_ops::
max_scale_trace(const vil_image_view<float> & input,float min_scale,float max_scale,float scale_inc)939 max_scale_trace(const vil_image_view<float>& input,
940 float min_scale, float max_scale, float scale_inc)
941 {
942 unsigned ni = input.ni(), nj = input.nj();
943 vil_image_view<float> tr_max, sc;
944 tr_max.set_size(ni, nj);
945 tr_max.fill(0.0f);
946 sc.set_size(ni, nj);
947 sc.fill(min_scale);
948 for (float s = min_scale; s<=max_scale; s+=scale_inc)
949 {
950 vil_image_view<float> smooth = brip_vil_float_ops::gaussian(input, s);
951 auto N = static_cast<unsigned>(2.0f*s);
952 vil_image_view<float> tr =
953 brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
954 for (unsigned r = 0; r<nj; ++r)
955 for (unsigned c = 0; c<ni; ++c)
956 {
957 float trv = s*s*tr(c,r);
958 if (trv>tr_max(c,r))
959 {
960 tr_max(c,r) = trv;
961 sc(c,r) = s;
962 }
963 }
964 }
965 return sc;
966 }
967
968 //: exactly same as max_scale_trace, only return the image with actual trace values at max scales instead of the image with max scale values
969 vil_image_view<float> brip_vil_float_ops::
max_scale_trace_value(const vil_image_view<float> & input,float min_scale,float max_scale,float scale_inc)970 max_scale_trace_value(const vil_image_view<float>& input,
971 float min_scale, float max_scale, float scale_inc)
972 {
973 unsigned ni = input.ni(), nj = input.nj();
974 vil_image_view<float> tr_max, sc;
975 tr_max.set_size(ni, nj);
976 tr_max.fill(0.0f);
977 sc.set_size(ni, nj);
978 sc.fill(min_scale);
979 for (float s = min_scale; s<=max_scale; s+=scale_inc)
980 {
981 vil_image_view<float> smooth = brip_vil_float_ops::gaussian(input, s);
982 auto N = static_cast<unsigned>(2.0f*s);
983 vil_image_view<float> tr =
984 brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
985 for (unsigned r = 0; r<nj; ++r)
986 for (unsigned c = 0; c<ni; ++c)
987 {
988 float trv = s*s*tr(c,r);
989 if (trv>tr_max(c,r))
990 {
991 tr_max(c,r) = trv;
992 sc(c,r) = s;
993 }
994 }
995 }
996 // mask the region where integration region extends outside the image borders
997 auto N = static_cast<unsigned>(5.0f*max_scale);
998 unsigned njmax = nj-N;
999 unsigned nimax = ni-N;
1000 for (unsigned r = 0; r<nj; ++r)
1001 for (unsigned c = 0; c<ni; ++c)
1002 {
1003 if (r <= N || r >= njmax || c <= N || c >= nimax)
1004 tr_max(c,r) = 0.0f;
1005 }
1006
1007 // normalize the trace values
1008 float min_b,max_b;
1009 vil_math_value_range(tr_max,min_b,max_b);
1010 std::cout << "in trace max image min value: " << min_b << " max value: " << max_b << std::endl;
1011
1012 vil_image_view<double> tr_normalized, tr_stretched;
1013 vil_convert_stretch_range(tr_max, tr_normalized, 0.0f, 1.0f);
1014 vil_convert_stretch_range(tr_max, tr_stretched, 0.0f, 256.0f);
1015 vil_image_view<float> tr_cast;
1016 vil_convert_cast(tr_stretched, tr_cast);
1017 vil_image_view<vxl_byte> tr_cast_byte;
1018 vil_convert_cast(tr_stretched, tr_cast_byte);
1019 vil_save_image_resource(vil_new_image_resource_of_view(tr_cast_byte), "D:\\projects\\vehicle_rec_on_models\\trace_image.png");
1020
1021 #if 0 // commented out
1022 // investigate illumination invariance
1023 vil_image_view<float> modified_input = input;
1024 vil_image_view<float> const_image(ni, nj);
1025 const_image.fill(3.0f);
1026 for (unsigned r = 0; r<nj; ++r)
1027 for (unsigned c = 0; c<ni; ++c)
1028 {
1029 if (modified_input(c,r) < 120 || modified_input(c,r) > 95) = modified_input(c,r)*const_image(c,r);
1030 }
1031
1032 vil_image_view<vxl_byte> mod_cast;
1033 vil_convert_cast(modified_input, mod_cast);
1034 vil_save_image_resource(vil_new_image_resource_of_view(mod_cast), "D:\\projects\\vehicle_rec_on_models\\modified_image.png");
1035
1036 vil_image_view<float> tr_max2, sc2;
1037 tr_max2.set_size(ni, nj);
1038 tr_max2.fill(0.0f);
1039 sc2.set_size(ni, nj);
1040 sc2.fill(min_scale);
1041 for (float s = min_scale; s<=max_scale; s+=scale_inc)
1042 {
1043 vil_image_view<float> smooth = brip_vil_float_ops::gaussian(modified_input, s);
1044 unsigned N = static_cast<unsigned>(2.0f*s);
1045 vil_image_view<float> tr =
1046 brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
1047 for (unsigned r = 0; r<nj; ++r)
1048 for (unsigned c = 0; c<ni; ++c)
1049 {
1050 float trv = s*s*tr(c,r);
1051 if (trv>tr_max2(c,r))
1052 {
1053 tr_max2(c,r) = trv;
1054 sc2(c,r) = s;
1055 }
1056 }
1057 }
1058 // mask the region where integration region extends outside the image borders
1059 for (unsigned r = 0; r<nj; ++r)
1060 for (unsigned c = 0; c<ni; ++c)
1061 {
1062 if (r <= N || r >= njmax || c <= N || c >= nimax)
1063 tr_max2(c,r) = 0.0f;
1064 }
1065
1066 vil_math_value_range(tr_max2,min_b,max_b);
1067 std::cout << "in trace max2 image min value: " << min_b << " max value: " << max_b << std::endl;
1068
1069 // normalize
1070 vil_image_view<double> tr_normalized2;
1071 vil_convert_stretch_range(tr_max2, tr_normalized2, 0.0f, 1.0f);
1072
1073 vil_image_view<double> difference_image = tr_normalized;
1074 for (unsigned r = 0; r<nj; ++r)
1075 for (unsigned c = 0; c<ni; ++c)
1076 {
1077 difference_image(c,r) = std::abs(difference_image(c,r) - tr_normalized2(c,r));
1078 }
1079 double min_bd, max_bd;
1080 vil_math_value_range(difference_image,min_bd,max_bd);
1081 std::cout << "in difference image of normalized trace images min value: " << min_bd << " max value: " << max_bd << std::endl;
1082 vil_image_view<vxl_byte> dif_cast;
1083 vil_convert_cast(difference_image, dif_cast);
1084 vil_save_image_resource(vil_new_image_resource_of_view(dif_cast), "D:\\projects\\vehicle_rec_on_models\\difference_image.png");
1085
1086 vil_image_view<double> tr_stretched2;
1087 vil_convert_stretch_range(tr_max2, tr_stretched2, 0.0f, 256.0f);
1088 vil_image_view<float> tr_cast2;
1089 vil_convert_cast(tr_stretched2, tr_cast2);
1090 vil_image_view<vxl_byte> tr_cast2_byte;
1091 vil_convert_cast(tr_stretched2, tr_cast2_byte);
1092 vil_save_image_resource(vil_new_image_resource_of_view(tr_cast2_byte), "D:\\projects\\vehicle_rec_on_models\\trace_image2.png");
1093
1094 return tr_cast2;
1095 #else // 0
1096 return tr_cast;
1097 #endif // 0
1098 }
1099
1100
1101 //---------------------------------------------------------------------
1102 // Lucas-Kanade motion vectors: Solve for the motion vectors over a
1103 // (2n+1)x(2n+1) neighborhood. The time derivative of intensity is computed
1104 // from the previous_frame. The threshold eliminates small values of
1105 // the product of the time derivative and the motion matrix eigenvalues,
1106 // i.e, |lambda_1*lambda_2*dI/dt|<thresh. Thus motion is only reported when
1107 // the solution is well-conditioned.
1108 //
1109 void
Lucas_KanadeMotion(vil_image_view<float> & current_frame,vil_image_view<float> & previous_frame,int n,double thresh,vil_image_view<float> & vx,vil_image_view<float> & vy)1110 brip_vil_float_ops::Lucas_KanadeMotion(vil_image_view<float> & current_frame,
1111 vil_image_view<float> & previous_frame,
1112 int n, double thresh,
1113 vil_image_view<float>& vx,
1114 vil_image_view<float>& vy)
1115 {
1116 int N = (2*n+1)*(2*n+1);
1117 int w = static_cast<int>(current_frame.ni()), h = static_cast<int>(current_frame.nj());
1118 vil_image_view<float> grad_x, grad_y, diff;
1119 grad_x.set_size(w,h);
1120 grad_y.set_size(w,h);
1121 // compute the gradient vector and the time derivative
1122 brip_vil_float_ops::gradient_3x3(current_frame, grad_x, grad_y);
1123 diff = brip_vil_float_ops::difference(previous_frame, current_frame);
1124 vul_timer t;
1125 // sum the motion terms over the (2n+1)x(2n+1) neighborhood.
1126 for (int y = n; y<h-n;y++)
1127 for (int x = n; x<w-n;x++)
1128 {
1129 float IxIx=0, IxIy=0, IyIy=0, IxIt=0, IyIt=0;
1130 for (int i = -n; i<=n; i++)
1131 for (int j = -n; j<=n; j++)
1132 {
1133 float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
1134 float dt = diff(x+i, y+j);
1135 IxIx += gx*gx;
1136 IxIy += gx*gy;
1137 IyIy += gy*gy;
1138 IxIt += gx*dt;
1139 IyIt += gy*dt;
1140 }
1141 //Divide by the number of pixels in the neighborhood
1142 IxIx/=N; IxIy/=N; IyIy/=N; IxIt/=N; IyIt/=N;
1143 auto det = float(IxIx*IyIy-IxIy*IxIy);
1144 //Eliminate small motion factors
1145 float dif = diff(x,y);
1146 float motion_factor = std::fabs(det*dif);
1147 if (motion_factor<thresh)
1148 {
1149 vx(x,y) = 0.0f;
1150 vy(x,y) = 0.0f;
1151 continue;
1152 }
1153 //solve for the motion vector
1154 vx(x,y) = (IyIy*IxIt-IxIy*IyIt)/det;
1155 vy(x,y) = (-IxIy*IxIt + IxIx*IyIt)/det;
1156 }
1157 brip_vil_float_ops::fill_x_border(vx, n, 0.0f);
1158 brip_vil_float_ops::fill_y_border(vx, n, 0.0f);
1159 brip_vil_float_ops::fill_x_border(vy, n, 0.0f);
1160 brip_vil_float_ops::fill_y_border(vy, n, 0.0f);
1161 #ifdef DEBUG
1162 std::cout << "\nCompute Lucas-Kanade in " << t.real() << " msecs.\n";
1163 #endif
1164 }
1165
1166 //: computes Lucas-Kanade optical flow on the complete input views
1167 void brip_vil_float_ops::
lucas_kanade_motion_on_view(vil_image_view<float> const & curr_frame,vil_image_view<float> const & prev_frame,const double thresh,float & vx,float & vy)1168 lucas_kanade_motion_on_view(vil_image_view<float> const& curr_frame,
1169 vil_image_view<float> const& prev_frame,
1170 const double thresh,
1171 float& vx,
1172 float& vy)
1173 {
1174 unsigned w = curr_frame.ni(), h = curr_frame.nj();
1175 unsigned N = w*h;
1176 vil_image_view<float> grad_x, grad_y;
1177 grad_x.set_size(w,h);
1178 grad_y.set_size(w,h);
1179 // compute the gradient vector and the time derivative
1180 brip_vil_float_ops::gradient_3x3(curr_frame, grad_x, grad_y);
1181 vil_image_view<float> diff =
1182 brip_vil_float_ops::difference(prev_frame, curr_frame);
1183
1184 // sum the motion terms over the view
1185 float IxIx=0.0f, IxIy=0.0f, IyIy=0.0f, IxIt=0.0f, IyIt=0.0f, dsum=0.0f;
1186 for (unsigned j = 0; j<h; j++)
1187 for (unsigned i = 0; i<w; i++)
1188 {
1189 float gx = grad_x(i, j), gy = grad_y(i, j);
1190 float dt = diff(i, j);
1191 dsum += dt*dt;
1192 IxIx += gx*gx;
1193 IxIy += gx*gy;
1194 IyIy += gy*gy;
1195 IxIt += gx*dt;
1196 IyIt += gy*dt;
1197 }
1198 // Divide by the number of pixels in the neighborhood
1199 IxIx/=(float)N; IxIy/=(float)N; IyIy/=(float)N; IxIt/=(float)N; IyIt/=(float)N; dsum/=(float)N;
1200 auto det = float(IxIx*IyIy-IxIy*IxIy);
1201 // Eliminate small motion factors
1202 float dif = std::sqrt(dsum);
1203 float motion_factor = std::fabs(det*dif);
1204 if (motion_factor<thresh)
1205 {
1206 vx = 0.0f;
1207 vy = 0.0f;
1208 return;
1209 }
1210 // solve for the motion vector
1211 vx = (IyIy*IxIt-IxIy*IyIt)/det;
1212 vy = (-IxIy*IxIt + IxIx*IyIt)/det;
1213 }
1214
1215 //------------------------------------------------------------------------
1216 // Assume that curr match region is larger than prev_region by the required
1217 // search ranges. Step through the search and output the shift that
1218 // maximizes the correlation. zero_i and zero_j indicate the curr_image
1219 // pixel location corresponding to no velocity between frames
1220 void brip_vil_float_ops::
velocity_by_correlation(vil_image_view<float> const & curr_image,vil_image_view<float> const & prev_region,const unsigned start_i,const unsigned end_i,const unsigned start_j,const unsigned end_j,const unsigned zero_i,const unsigned zero_j,float & vx,float & vy)1221 velocity_by_correlation(vil_image_view<float> const& curr_image,
1222 vil_image_view<float> const& prev_region,
1223 const unsigned start_i, const unsigned end_i,
1224 const unsigned start_j, const unsigned end_j,
1225 const unsigned zero_i, const unsigned zero_j,
1226 float& vx,
1227 float& vy)
1228 {
1229 unsigned ni = prev_region.ni(), nj = prev_region.nj();
1230 float corr_max = -10.0f;
1231 auto vx0 = static_cast<float>(zero_i);
1232 auto vy0 = static_cast<float>(zero_j);
1233 auto area = static_cast<float>(ni*nj);
1234 vx = 0; vy = 0;
1235 unsigned max_i = start_i, max_j = start_j;
1236 for (unsigned j = start_j; j<=end_j; ++j)
1237 for (unsigned i = start_i; i<=end_i; ++i)
1238 {
1239 float si1 = 0, si2 = 0, si1i1 = 0, si2i2 = 0, si1i2 = 0;
1240 //sum over the region
1241 for (unsigned r = 0; r<nj; ++r)
1242 for (unsigned c = 0; c<ni; ++c)
1243 {
1244 float I1 = prev_region(c, r);
1245 float I2 = curr_image(i+c, j+r);
1246 si1 += I1; si2 += I2;
1247 si1i1 += I1*I1;
1248 si2i2 += I2*I2;
1249 si1i2 += I1*I2;
1250 }
1251 float corr = cross_corr(area, si1, si2, si1i1, si2i2,si1i2, 1.0f);
1252 if (corr>corr_max)
1253 {
1254 corr_max = corr;
1255 max_i = i; max_j = j;
1256 }
1257 #if 0
1258 float di = i-vx0, dj = j-vy0;
1259 std::cout << di << '\t' << dj << '\t' << corr << '\n';
1260 #endif
1261 }
1262 // the velocity is given by the max indices relative to the zero location
1263 vx = static_cast<float>(max_i)- vx0;
1264 vy = static_cast<float>(max_j) - vy0;
1265 // std::cout << '(' << vx << ' ' << vy << "): " << corr_max << '\n';
1266 }
1267
1268 //---------------------------------------------------------------------
1269 // Horn-Schunk method for calc. motion vectors: Solve for the motion vectors
1270 // iteratively using a cost function with two terms (RHS of optical flow eqn
1271 // and the magnitude of spatial derivatives of the velocity field
1272 // (so that pixel-to-pixel variations are small). The second term is
1273 // weighted by alpha_coef term. The iteration goes on until the error
1274 // reaches below err_thresh
1275 // Error conditions:
1276 // - -1 - current_frame and previous_frame are equal
1277 // - -2 - at least one input frame or internal process image is all zeros
1278 // - 0 - routine was successful
1279
1280 int brip_vil_float_ops::
Horn_SchunckMotion(vil_image_view<float> const & current_frame,vil_image_view<float> const & previous_frame,vil_image_view<float> & vx,vil_image_view<float> & vy,const float alpha_coef,const int no_of_iterations)1281 Horn_SchunckMotion(vil_image_view<float> const& current_frame,
1282 vil_image_view<float> const& previous_frame,
1283 vil_image_view<float>& vx,
1284 vil_image_view<float>& vy,
1285 const float alpha_coef,
1286 const int no_of_iterations)
1287 {
1288 // Check for equal images
1289 if (vil_image_view_deep_equality (previous_frame, current_frame ) )
1290 {
1291 std::cout<<"Images are same";
1292 return -1;
1293 }
1294
1295 // Declarations
1296 unsigned w = current_frame.ni(), h = current_frame.nj();
1297
1298 vil_image_view<float> grad_x, grad_y, diff;
1299
1300 vil_image_view<float> temp1;
1301 vil_image_view<float> temp2;
1302
1303 vil_image_view<float> emptyimg;
1304
1305 // Size Init
1306
1307 grad_x.set_size(w,h);
1308 grad_y.set_size(w,h);
1309 diff.set_size(w,h);
1310 temp1.set_size(w,h);
1311 temp2.set_size(w,h);
1312
1313 emptyimg.set_size(w,h);
1314
1315 temp1.fill(0.0);
1316 temp2.fill(0.0);
1317
1318 // Initialization
1319 for (unsigned y = 0; y<h; y++)
1320 for (unsigned x = 0; x<w; x++)
1321 {
1322 vx(x,y)=0.0f;
1323 vy(x,y)=0.0f;
1324 diff (x,y)=0.0f;
1325 grad_x (x, y)= 0.0f;
1326 grad_y (x, y)= 0.0f;
1327
1328 emptyimg (x, y) = 0.0f;
1329 }
1330
1331 // Check for empty images
1332 if ( (vil_image_view_deep_equality (emptyimg, current_frame )) || (vil_image_view_deep_equality(emptyimg, previous_frame)))
1333 {
1334 std::cout<<"Image is empty";
1335 return -2;
1336 }
1337
1338 // compute the gradient vector for current and previous
1339 brip_vil_float_ops::gradient_3x3 (current_frame , grad_x , grad_y);
1340 brip_vil_float_ops::gradient_3x3 (previous_frame , temp1 , temp2);
1341
1342 // Grad = 0.5* Grad(current) + 0.5 * Grad(previous)
1343 vil_math_add_image_fraction(grad_x, 0.5, temp1, 0.5);
1344 vil_math_add_image_fraction(grad_y, 0.5, temp2, 0.5);
1345 if ( (vil_image_view_deep_equality(emptyimg, grad_x)) ||
1346 (vil_image_view_deep_equality(emptyimg, grad_y)) )
1347 {
1348 std::cout<<"Gradient Image is empty";
1349 return -2;
1350 }
1351
1352 temp1.fill(0.0);
1353 temp2.fill(0.0);
1354
1355 // Average the local intensities over 3x3 region
1356 temp1 = brip_vil_float_ops::average_NxN(previous_frame, 3);
1357 temp2 = brip_vil_float_ops::average_NxN(current_frame, 3);
1358 if (vil_image_view_deep_equality(emptyimg, temp1) ||
1359 vil_image_view_deep_equality(emptyimg, temp2))
1360 {
1361 std::cout<<"Averaged Image is empty";
1362 return -2;
1363 }
1364
1365 // Compute the time derivative (difference of local average intensities)
1366 // diff = dI/dt
1367 diff = brip_vil_float_ops::difference(temp1 , temp2);
1368 if (vil_image_view_deep_equality(emptyimg, diff) )
1369 {
1370 std::cout<<"Difference Image is empty";
1371 return -2;
1372 }
1373
1374 temp1.fill(0.0);
1375 temp2.fill(0.0);
1376 // Iterate
1377 #ifdef DEBUG
1378 vul_timer t;
1379 #endif
1380 for (int i=0;i<no_of_iterations;i++)
1381 {
1382 // Update vx and vy
1383 //Smoothed velocities on 3x3 region
1384 temp1 = brip_vil_float_ops::average_NxN (vx, 3);
1385 temp2 = brip_vil_float_ops::average_NxN (vy, 3);
1386 for (unsigned y = 1; y+1<h; ++y)
1387 for (unsigned x = 1; x+1<w; ++x)
1388 {
1389 float tempx = temp1(x,y);
1390 float tempy = temp2(x,y);
1391
1392 float gx = grad_x(x, y), gy = grad_y(x, y);
1393
1394 float dt = diff(x, y);
1395 // _____
1396 // term = (v(x,y).Grad(x,y) + dI/dt(x,y))/(alpha + |Grad(x,y)|^2)
1397 // term is the brightness constraint normalized by gradient mag.
1398 //
1399 float term =
1400 ( (gx * tempx) + (gy * tempy) + dt )/ (alpha_coef + gx*gx + gy*gy);
1401
1402 // ______
1403 //v(x,y) = v(x,y) - Grad(x,y)* term
1404 vx(x,y) = tempx - (gx * term);
1405 vy(x,y) = tempy - (gy * term);
1406 }
1407
1408 #ifdef DEBUG
1409 std::cout << "Iteration No " << i << '\n';
1410 #endif
1411 brip_vil_float_ops::fill_x_border(vx, 1, 0.0f);
1412 brip_vil_float_ops::fill_y_border(vx, 1, 0.0f);
1413 brip_vil_float_ops::fill_x_border(vy, 1, 0.0f);
1414 brip_vil_float_ops::fill_y_border(vy, 1, 0.0f);
1415 }
1416 #ifdef DEBUG
1417 std::cout << "\nCompute Horn-Schunck iteration in " << t.real() << " msecs.\n";
1418 #endif
1419 return 0;
1420 }
1421
fill_x_border(vil_image_view<float> & image,unsigned w,float value)1422 void brip_vil_float_ops::fill_x_border(vil_image_view<float> & image,
1423 unsigned w, float value)
1424 {
1425 unsigned width = image.ni(), height = image.nj();
1426 if (2*w>width)
1427 {
1428 std::cout << "In brip_vil_float_ops::fill_x_border(..) - 2xborder exceeds image width\n";
1429 return;
1430 }
1431 for (unsigned y = 0; y<height; y++)
1432 for (unsigned x = 0; x<w; x++)
1433 image(x, y) = value;
1434
1435 for (unsigned y = 0; y<height; y++)
1436 for (unsigned x = width-w; x<width; x++)
1437 image(x, y) = value;
1438 }
1439
fill_y_border(vil_image_view<float> & image,unsigned h,float value)1440 void brip_vil_float_ops::fill_y_border(vil_image_view<float> & image,
1441 unsigned h, float value)
1442 {
1443 unsigned width = image.ni(), height = image.nj();
1444 if (2*h>height)
1445 {
1446 std::cout << "In brip_vil_float_ops::fill_y_border(..) - 2xborder exceeds image height\n";
1447 return;
1448 }
1449 for (unsigned y = 0; y<h; y++)
1450 for (unsigned x = 0; x<width; x++)
1451 image(x, y) = value;
1452
1453 for (unsigned y = height-h; y<height; y++)
1454 for (unsigned x = 0; x<width; x++)
1455 image(x, y) = value;
1456 }
1457
1458 vil_image_view<vxl_byte>
convert_to_byte(vil_image_view<float> const & image)1459 brip_vil_float_ops::convert_to_byte(vil_image_view<float> const& image)
1460 {
1461 // determine the max min values
1462 float min_val = vnl_numeric_traits<float>::maxval;
1463 float max_val = -min_val;
1464 unsigned w = image.ni(), h = image.nj();
1465 vil_image_view<vxl_byte> output;
1466 output.set_size(w,h);
1467 for (unsigned y = 0; y<h; y++)
1468 for (unsigned x = 0; x<w; x++)
1469 {
1470 min_val = std::min(min_val, image(x,y));
1471 max_val = std::max(max_val, image(x,y));
1472 }
1473 float range = max_val-min_val;
1474 if (range == 0.f)
1475 range = 1.f;
1476 else
1477 range = 255.f/range;
1478 for (unsigned y = 0; y<h; y++)
1479 for (unsigned x = 0; x<w; x++)
1480 {
1481 float v = (image(x,y)-min_val)*range;
1482 if(v>255.0f) v = 255.0f;//could have small round off above 255
1483 output(x,y) = (unsigned char)v;
1484 }
1485 return output;
1486 }
1487
1488 //------------------------------------------------------------
1489 //: Convert the range between min_val and max_val to 255
1490 vil_image_view<vxl_byte>
convert_to_byte(vil_image_view<float> const & image,float min_val,float max_val)1491 brip_vil_float_ops::convert_to_byte(vil_image_view<float> const& image,
1492 float min_val, float max_val)
1493 {
1494 unsigned w = image.ni(), h = image.nj();
1495 vil_image_view<vxl_byte> output;
1496 output.set_size(w,h);
1497 float range = max_val-min_val;
1498 if (range == 0.f)
1499 range = 1.f;
1500 else
1501 range = 256.f/range;
1502 for (unsigned y = 0; y<h; y++)
1503 for (unsigned x = 0; x<w; x++)
1504 {
1505 float v = (image(x,y)-min_val)*range;
1506 if (v>=256)
1507 v=255;
1508 else if (v<0)
1509 v=0;
1510 output(x,y) = (unsigned char)v;
1511 }
1512
1513 return output;
1514 }
1515
1516 vil_image_view<vxl_byte> brip_vil_float_ops::
convert_to_byte(vil_image_view<vxl_uint_16> const & image,vxl_uint_16 min_val,vxl_uint_16 max_val)1517 convert_to_byte(vil_image_view<vxl_uint_16> const& image,
1518 vxl_uint_16 min_val, vxl_uint_16 max_val)
1519 {
1520 unsigned ni = image.ni(), nj = image.nj();
1521 vil_image_view<vxl_byte> output;
1522 output.set_size(ni, nj);
1523 auto range = static_cast<float>(max_val-min_val);
1524 if (!range)
1525 range = 1.f;
1526 else
1527 range = 256.f/range;
1528 for (unsigned r = 0; r<nj; r++)
1529 for (unsigned c = 0; c<ni; c++)
1530 {
1531 float v = float(image(c, r)-min_val)*range;
1532 if (v>=256)
1533 v=255;
1534 else if (v<0)
1535 v=0;
1536 output(c, r) = static_cast<unsigned char>(v);
1537 }
1538 return output;
1539 }
1540
1541 // Note this is a more standard interface than convert_to_grey
1542 vil_image_view<vxl_byte>
convert_to_byte(vil_image_resource_sptr const & image)1543 brip_vil_float_ops::convert_to_byte(vil_image_resource_sptr const& image)
1544 {
1545 return brip_vil_float_ops::convert_to_grey(*image);
1546 }
1547
1548 vil_image_view<vxl_uint_16>
convert_to_short(vil_image_view<float> const & image,float min_val,float max_val)1549 brip_vil_float_ops::convert_to_short(vil_image_view<float> const& image,
1550 float min_val, float max_val)
1551 {
1552 unsigned w = image.ni(), h = image.nj();
1553 float max_short = 65355.f;
1554 vil_image_view<vxl_uint_16> output;
1555 output.set_size(w,h);
1556 float range = max_val-min_val;
1557 if (!range)
1558 range = 1;
1559 else
1560 range = max_short/range;
1561 for (unsigned y = 0; y<h; y++)
1562 for (unsigned x = 0; x<w; x++)
1563 {
1564 float v = (image(x,y)-min_val)*range;
1565 if (v>max_short)
1566 v=max_short;
1567 if (v<0)
1568 v=0;
1569 output(x,y) = (unsigned short)v;
1570 }
1571 return output;
1572 }
1573
1574 //: converts a float image to an unsigned short (16-bit) image.
1575 // range determined automatically
1576 vil_image_view<vxl_uint_16>
convert_to_short(vil_image_view<float> const & image)1577 brip_vil_float_ops::convert_to_short(vil_image_view<float> const& image)
1578 {
1579 float minv = vnl_numeric_traits<float>::maxval;
1580 float maxv = -minv;
1581 unsigned ni = image.ni(), nj = image.nj();
1582 for (unsigned j = 0; j<nj; ++j)
1583 for (unsigned i = 0; i<ni; ++i)
1584 {
1585 float v = image(i,j);
1586 if (v<minv)
1587 minv = v;
1588 if (v>maxv)
1589 maxv = v;
1590 }
1591 return brip_vil_float_ops::convert_to_short(image, minv, maxv);
1592 }
1593
1594 vil_image_view<vxl_uint_16>
convert_to_short(vil_image_resource_sptr const & image)1595 brip_vil_float_ops::convert_to_short(vil_image_resource_sptr const& image)
1596 {
1597 // Check if the image is a float
1598 if (image->nplanes()==1 &&image->pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
1599 {
1600 vil_image_view<float> temp = image->get_view();
1601 float vmin=0, vmax= 65355;
1602 vil_math_value_range<float>(temp, vmin, vmax);
1603 return brip_vil_float_ops::convert_to_short(temp, vmin, vmax);
1604 }
1605
1606 // Here we assume that the image is an unsigned char
1607 if (image->nplanes()==1&&image->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
1608 {
1609 vil_image_view<unsigned char > temp = image->get_view();
1610 vil_image_view<unsigned short> short_image;
1611 unsigned width = temp.ni(), height = temp.nj();
1612 short_image.set_size(width, height);
1613 for (unsigned y = 0; y<height; y++)
1614 for (unsigned x = 0; x<width; x++)
1615 short_image(x,y) = static_cast<unsigned short>(temp(x,y));
1616 return temp;
1617 }
1618
1619 // Here the image is an unsigned short (16-bit) image so just return it
1620 if (image->nplanes()==1&&image->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
1621 {
1622 vil_image_view<unsigned short > temp = image->get_view();
1623 return temp;
1624 }
1625
1626 // the image is color so we should convert it to greyscale
1627 // Here we assume the color elements are unsigned char.
1628 if (image->nplanes()==3&&image->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
1629 {
1630 vil_image_view<vxl_byte> color_image = image->get_view();
1631 unsigned width = color_image.ni(), height = color_image.nj();
1632 // the output image
1633 vil_image_view<unsigned short> short_image;
1634 short_image.set_size(width, height);
1635 for (unsigned y = 0; y<height; y++)
1636 for (unsigned x = 0; x<width; x++)
1637 {
1638 double v = color_image(x,y,0)+color_image(x,y,1)+color_image(x,y,2);
1639 v/=3.0;
1640 short_image(x,y) = static_cast<unsigned short>(v);
1641 }
1642 return short_image;
1643 }
1644
1645 // the image is multispectral so we should convert it to greyscale
1646 // Here we assume the color elements are unsigned short (16-bit).
1647 if (image->nplanes()==4&&image->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
1648 {
1649 vil_image_view<unsigned short > mband_image = image->get_view();
1650 unsigned width = mband_image.ni(), height = mband_image.nj();
1651 // the output image
1652 vil_image_view<unsigned short> short_image;
1653 short_image.set_size(width, height);
1654 for (unsigned y = 0; y<height; y++)
1655 for (unsigned x = 0; x<width; x++)
1656 {
1657 unsigned short v = 0;
1658 for (unsigned p = 0; p<4; ++p)
1659 v += mband_image(x, y, p);
1660 v/=4;
1661 short_image(x,y) = v;
1662 }
1663 return short_image;
1664 }
1665
1666 // If we get here then the input is not a type we handle so return a null view
1667 return vil_image_view<vxl_uint_16>();
1668 }
1669
1670 vil_image_view<float>
convert_to_float(vil_image_view<vxl_byte> const & image)1671 brip_vil_float_ops::convert_to_float(vil_image_view<vxl_byte> const& image)
1672 {
1673 vil_image_view<float> output;
1674 unsigned ni = image.ni(), nj = image.nj(), np = image.nplanes();
1675 output.set_size(ni,nj, np);
1676 for (unsigned j = 0; j<nj; j++)
1677 for (unsigned i = 0; i<ni; i++)
1678 for (unsigned p = 0; p<np; ++p)
1679 output(i,j,p) = (float)image(i,j,p);
1680 return output;
1681 }
1682
1683 vil_image_view<float>
convert_to_float(vil_image_view<vxl_uint_16> const & image)1684 brip_vil_float_ops::convert_to_float(vil_image_view<vxl_uint_16> const& image)
1685 {
1686 vil_image_view<float> output;
1687 unsigned w = image.ni(), h = image.nj();
1688 output.set_size(w,h);
1689 for (unsigned y = 0; y<h; y++)
1690 for (unsigned x = 0; x<w; x++)
1691 output(x,y) = (float)image(x,y);
1692 return output;
1693 }
1694
1695 vil_image_view<bool>
convert_to_bool(vil_image_view<vxl_byte> const & image)1696 brip_vil_float_ops::convert_to_bool(vil_image_view<vxl_byte> const& image)
1697 {
1698 vil_image_view<bool> output;
1699 unsigned w = image.ni(), h = image.nj();
1700 output.set_size(w,h);
1701 for (unsigned y = 0; y<h; y++)
1702 for (unsigned x = 0; x<w; x++)
1703 if (image(x,y) > 128)
1704 output(x,y)=true;
1705 else
1706 output(x,y)=false;
1707 return output;
1708 }
1709
1710 vil_image_view<float>
convert_to_float(vil_image_view<vil_rgb<vxl_byte>> const & image)1711 brip_vil_float_ops::convert_to_float(vil_image_view<vil_rgb<vxl_byte> > const& image)
1712 {
1713 vil_image_view<float> output;
1714 unsigned w = image.ni(), h = image.nj();
1715 output.set_size(w,h);
1716 for (unsigned y = 0; y<h; y++)
1717 for (unsigned x = 0; x<w; x++)
1718 {
1719 vil_rgb<vxl_byte> rgb = image(x,y);
1720 output(x,y) = (float)rgb.grey();
1721 }
1722 return output;
1723 }
1724
rgb_to_ihs(vil_rgb<vxl_byte> const & rgb,float & i,float & h,float & s)1725 void brip_vil_float_ops::rgb_to_ihs(vil_rgb<vxl_byte> const& rgb,
1726 float& i, float& h, float& s)
1727 {
1728 float r=rgb.R();
1729 float g=rgb.G();
1730 float b=rgb.B();
1731
1732 float maxval = std::max(r,std::max(g,b));
1733 float minval = std::min(r,std::min(g,b));
1734
1735 float delta = maxval - minval;
1736 i = maxval;
1737 if (maxval == 0)
1738 s = 0;
1739 else
1740 s = delta / maxval;
1741
1742 if (s==0)
1743 h = 0; //!< (Hue is undefined)
1744
1745 if (r== maxval)
1746 h = (g - b) / delta ; //!< (between yellow and magenta)
1747 if (g == maxval)
1748 h = 2 + (b - r)/delta ; //!< (between cyan and yellow)
1749 if (b == maxval)
1750 h = 4 + (r - g) / delta; //!< (between magenta and cyan)
1751 h *= 60.f; //!< (convert Hue to degrees)
1752 if (h < 0.f)
1753 h += 360.f; //!< (Hue must be positive)
1754 else if (h >= 360.f)
1755 h -= 360.f; //!< (Hue must be less than 360)
1756
1757 h *= 256.0f / 360.0f; // range 0 .. 255
1758 s *= 256.0f; // range 0 .. 255
1759 }
1760
rgb_to_ihs_tsai(vil_rgb<vxl_byte> const & rgb,float & i,float & h,float & s)1761 void brip_vil_float_ops::rgb_to_ihs_tsai(vil_rgb<vxl_byte> const& rgb, float& i, float& h, float& s)
1762 {
1763 float r=rgb.R();
1764 float g=rgb.G();
1765 float b=rgb.B();
1766 float sq6 = std::sqrt(6.0f);
1767 float V1, V2;
1768 i = (r+g+b)/3.0f;
1769 V1 = (2.0f*b-r-g)/sq6;
1770 V2 = (r-2.0f*g)/sq6;
1771 s = std::sqrt(V1*V1 + V2*V2);
1772 auto two_pi = static_cast<float>(vnl_math::twopi);
1773 float a = std::atan2(V2, V1)/two_pi;
1774 if (a<0.0f) a += 1.0f; // range [0..1)
1775 h = a*256.0f; // range [0..256)
1776 }
1777
ihs_to_rgb(vil_rgb<vxl_byte> & rgb,const float i,const float h,const float s)1778 void brip_vil_float_ops::ihs_to_rgb(vil_rgb<vxl_byte> & rgb,
1779 const float i, const float h, const float s)
1780 {
1781 // Reference: page 593 of Foley & van Dam
1782 float R = 0.0f;
1783 float G = 0.0f;
1784 float B = 0.0f;
1785
1786 if (s == 0) {
1787 R=i;
1788 G=i;
1789 B=i;
1790 }
1791 else if (s > 0.0)
1792 {
1793 float ss = s, hh = h;
1794 ss *= 1.f / 256.f;
1795 hh *= 6.f / 256.f;
1796
1797 float J = std::floor(hh);
1798 float F = hh - J;
1799 float P =( i * (1 - ss));
1800 float Q = (i * (1 - (ss * F)));
1801 float T = (i * (1 - (ss * (1 - F))));
1802
1803 if (J == 0) { R=i; G=T; B=P; }
1804 if (J == 1) { R=Q; G=i; B=P; }
1805 if (J == 2) { R=P; G=i; B=T; }
1806 if (J == 3) { R=P; G=Q; B=i; }
1807 if (J == 4) { R=T; G=P; B=i; }
1808 if (J == 5) { R=i; G=P; B=Q; }
1809 }
1810 rgb.r = (vxl_byte)R;
1811 rgb.g = (vxl_byte)G;
1812 rgb.b = (vxl_byte)B;
1813 }
1814
1815 void brip_vil_float_ops::
convert_to_IHS(vil_image_view<vil_rgb<vxl_byte>> const & image,vil_image_view<float> & I,vil_image_view<float> & H,vil_image_view<float> & S)1816 convert_to_IHS(vil_image_view<vil_rgb<vxl_byte> > const& image,
1817 vil_image_view<float>& I,
1818 vil_image_view<float>& H,
1819 vil_image_view<float>& S)
1820 {
1821 unsigned w = image.ni(), h = image.nj();
1822 I.set_size(w,h);
1823 H.set_size(w,h);
1824 S.set_size(w,h);
1825 for (unsigned r = 0; r < h; r++)
1826 for (unsigned c = 0; c < w; c++)
1827 {
1828 float in, hue, sat;
1829 rgb_to_ihs(image(c,r), in, hue, sat);
1830 I(c,r) = in;
1831 H(c,r) = hue;
1832 S(c,r) = sat;
1833 }
1834 }
1835
1836 void brip_vil_float_ops::
convert_to_IHS(vil_image_view<vxl_byte> const & image,vil_image_view<float> & I,vil_image_view<float> & H,vil_image_view<float> & S)1837 convert_to_IHS(vil_image_view<vxl_byte> const& image,
1838 vil_image_view<float>& I,
1839 vil_image_view<float>& H,
1840 vil_image_view<float>& S)
1841 {
1842 unsigned w = image.ni(), h = image.nj();
1843 I.set_size(w,h);
1844 H.set_size(w,h);
1845 S.set_size(w,h);
1846 for (unsigned r = 0; r < h; r++)
1847 for (unsigned c = 0; c < w; c++)
1848 {
1849 float in, hue, sat;
1850 vil_rgb<vxl_byte> imint(image(c,r,0),image(c,r,1),image(c,r,2));
1851 rgb_to_ihs(imint, in, hue, sat);
1852 I(c,r) = in;
1853 H(c,r) = hue;
1854 S(c,r) = sat;
1855 }
1856 }
1857
1858 void brip_vil_float_ops::
convert_to_IHS_tsai(vil_image_view<vxl_byte> const & image,vil_image_view<float> & I,vil_image_view<float> & H,vil_image_view<float> & S,vil_image_view<float> & ratio)1859 convert_to_IHS_tsai(vil_image_view<vxl_byte> const& image,
1860 vil_image_view<float>& I,
1861 vil_image_view<float>& H,
1862 vil_image_view<float>& S,
1863 vil_image_view<float>& ratio)
1864 {
1865 unsigned w = image.ni(), h = image.nj();
1866 I.set_size(w,h);
1867 H.set_size(w,h);
1868 S.set_size(w,h);
1869 ratio.set_size(w,h);
1870 for (unsigned r = 0; r < h; r++)
1871 for (unsigned c = 0; c < w; c++)
1872 {
1873 float in, hue, sat;
1874 vil_rgb<vxl_byte> imint(image(c,r,0),image(c,r,1),image(c,r,2));
1875 rgb_to_ihs_tsai(imint, in, hue, sat);
1876 I(c,r) = in;
1877 H(c,r) = hue;
1878 S(c,r) = sat;
1879 ratio(c,r) = (hue+1.0f)/(in+1.0f);
1880 }
1881 }
1882
1883 #if 0 // commented out
1884 void brip_vil_float_ops::
1885 display_IHS_as_RGB(vil_image_view<float> const& I,
1886 vil_image_view<float> const& H,
1887 vil_image_view<float> const& S,
1888 vil_image_view<vil_rgb<vxl_byte> >& image)
1889 {
1890 unsigned w = I.ni(), h = I.nj();
1891 image.set_size(w,h);
1892 float s = 256.0f/360.0f;
1893 for (unsigned r = 0; r < h; r++)
1894 for (unsigned c = 0; c < w; c++)
1895 {
1896 float in = I(c,r);
1897 float hue = s * H(c,r);
1898 float sat = S(c,r);
1899 if (in<0) in = 0;
1900 if (sat<0) sat = 0;
1901 if (hue<0) hue = 0;
1902 if (in>=256) in = 255;
1903 if (hue>=256) hue = 255;
1904 if (sat>=256) sat = 255;
1905 image(c,r).r = (vxl_byte)in;
1906 image(c,r).g = (vxl_byte)hue;
1907 image(c,r).b = (vxl_byte)sat;
1908 }
1909 }
1910 #endif // 0
1911
1912 // map so that intensity is proportional to saturation and hue is color
1913 void brip_vil_float_ops::
display_IHS_as_RGB(vil_image_view<float> const & I,vil_image_view<float> const & H,vil_image_view<float> const & S,vil_image_view<vil_rgb<vxl_byte>> & image)1914 display_IHS_as_RGB(vil_image_view<float> const& I,
1915 vil_image_view<float> const& H,
1916 vil_image_view<float> const& S,
1917 vil_image_view<vil_rgb<vxl_byte> >& image)
1918 {
1919 unsigned w = I.ni(), h = I.nj();
1920 image.set_size(w,h);
1921
1922 const auto deg_to_rad = float(vnl_math::pi_over_180);
1923 for (unsigned r = 0; r < h; r++)
1924 for (unsigned c = 0; c < w; c++)
1925 {
1926 float hue = H(c,r);
1927 float sat = 2.f*S(c,r);
1928 if (sat<0)
1929 sat = 0.f;
1930 else if (sat>=256)
1931 sat = 255.999f;
1932 float ang = deg_to_rad*hue;
1933 float cs = std::cos(ang), si = std::fabs(std::sin(ang));
1934 float red=0.0f, blue=0.0f;
1935 float green = si*sat;
1936 if (cs>=0)
1937 red = cs*sat;
1938 else
1939 blue = sat*(-cs);
1940 image(c,r).r = (vxl_byte)red;
1941 image(c,r).g = (vxl_byte)green;
1942 image(c,r).b = (vxl_byte)blue;
1943 }
1944 }
1945
1946 vil_image_view<vil_rgb<vxl_byte> > brip_vil_float_ops::
combine_color_planes(vil_image_view<vxl_byte> const & R,vil_image_view<vxl_byte> const & G,vil_image_view<vxl_byte> const & B)1947 combine_color_planes(vil_image_view<vxl_byte> const& R,
1948 vil_image_view<vxl_byte> const& G,
1949 vil_image_view<vxl_byte> const& B)
1950 {
1951 unsigned w = R.ni(), h = R.nj();
1952 vil_image_view<vil_rgb<vxl_byte> > image(w,h);
1953 for (unsigned r = 0; r < h; r++)
1954 for (unsigned c = 0; c < w; c++)
1955 {
1956 image(c,r).r = R(c,r);
1957 image(c,r).g = G(c,r);
1958 image(c,r).b = B(c,r);
1959 }
1960 return image;
1961 }
1962
1963 vil_image_view<vil_rgb<vxl_byte> >
combine_color_planes(vil_image_resource_sptr const & R,vil_image_resource_sptr const & G,vil_image_resource_sptr const & B)1964 brip_vil_float_ops::combine_color_planes(vil_image_resource_sptr const& R,
1965 vil_image_resource_sptr const& G,
1966 vil_image_resource_sptr const& B)
1967 {
1968 vil_image_view<vil_rgb<vxl_byte> > view(0,0);
1969 if (!R||!G||!B)
1970 return view; // return an empty view
1971 // determine the union of all the resources
1972 vbl_bounding_box<unsigned int, 2> b;
1973 unsigned int r_ni = R->ni(), r_nj = R->nj();
1974 unsigned int g_ni = G->ni(), g_nj = G->nj();
1975 unsigned int b_ni = B->ni(), b_nj = B->nj();
1976 b.update(r_ni, r_nj);
1977 b.update(g_ni, g_nj);
1978 b.update(b_ni, b_nj);
1979 unsigned int n_i = b.xmax(), n_j = b.ymax();
1980 view.set_size(n_i, n_j);
1981 vil_rgb<vxl_byte> zero(0, 0, 0);
1982 vil_image_view<float> fR = brip_vil_float_ops::convert_to_float(R);
1983 vil_image_view<vxl_byte> cR = brip_vil_float_ops::convert_to_byte(fR);
1984 vil_image_view<float> fG = brip_vil_float_ops::convert_to_float(G);
1985 vil_image_view<vxl_byte> cG = brip_vil_float_ops::convert_to_byte(fG);
1986 vil_image_view<float> fB = brip_vil_float_ops::convert_to_float(B);
1987 vil_image_view<vxl_byte> cB = brip_vil_float_ops::convert_to_byte(fB);
1988 for (vxl_uint_16 j = 0; j<n_j; ++j)
1989 for (vxl_uint_16 i = 0; i<n_i; ++i)
1990 {
1991 vil_rgb<vxl_byte> v = zero;
1992 if (i<r_ni&&j<r_nj)
1993 v.r = cR(i,j);
1994 if (i<g_ni&&j<g_nj)
1995 v.g= cG(i,j);
1996 if (i<b_ni&&j<b_nj)
1997 v.b= cB(i,j);
1998 view(i,j)=v;
1999 }
2000 return view;
2001 }
2002
2003 vil_image_view<float>
convert_to_float(vil_image_resource const & image)2004 brip_vil_float_ops::convert_to_float(vil_image_resource const& image)
2005 {
2006 vil_image_view<float> fimg;
2007 vil_pixel_format fmt = image.pixel_format();
2008 if (vil_pixel_format_num_components(image.pixel_format())==1)
2009 {
2010 if (fmt==VIL_PIXEL_FORMAT_UINT_16)
2011 {
2012 vil_image_view<unsigned short> temp=image.get_view();
2013 fimg = brip_vil_float_ops::convert_to_float(temp);
2014 }
2015 else if (fmt==VIL_PIXEL_FORMAT_BYTE)
2016 {
2017 vil_image_view<unsigned char> temp=image.get_view();
2018 fimg = brip_vil_float_ops::convert_to_float(temp);
2019 }
2020 else if (fmt==VIL_PIXEL_FORMAT_FLOAT)
2021 return image.get_view();
2022 }
2023 else if (vil_pixel_format_num_components(fmt)==3)
2024 {
2025 vil_image_view<vil_rgb<vxl_byte> > temp= image.get_view();
2026 fimg = brip_vil_float_ops::convert_to_float(temp);
2027 }
2028 else
2029 {
2030 std::cout << "In brip_vil_float_ops::convert_to_float - input not color or grey\n";
2031 return vil_image_view<float>();
2032 }
2033 return fimg;
2034 }
2035
2036 //-----------------------------------------------------------------
2037 //: Convert any image to an unsigned_char image
2038 vil_image_view<vxl_byte>
convert_to_grey(vil_image_resource const & image)2039 brip_vil_float_ops::convert_to_grey(vil_image_resource const& image)
2040 {
2041 // Check if the image is a float
2042 if (image.nplanes()==1 &&image.pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
2043 {
2044 vil_image_view<float> temp = image.get_view();
2045 float vmin=0.f, vmax=256.f;
2046 vil_math_value_range<float>(temp, vmin, vmax);
2047 return brip_vil_float_ops::convert_to_byte(temp, vmin, vmax);
2048 }
2049 if (image.nplanes()==1 &&image.pixel_format()==VIL_PIXEL_FORMAT_BOOL)
2050 {
2051 vil_image_view<bool> temp = image.get_view();
2052 unsigned nj = temp.nj(), ni = temp.ni();
2053 vil_image_view<unsigned char> out(ni, nj);
2054 out.fill(0);
2055 for (unsigned j = 0; j<nj; ++j)
2056 for (unsigned i = 0; i<ni; ++i)
2057 if (temp(i,j)) out(i,j) = 255;
2058 return out;
2059 }
2060
2061 // Here we assume that the image is an unsigned char
2062 // In this case we should just return it.
2063 if (image.nplanes()==1&&image.pixel_format()==VIL_PIXEL_FORMAT_BYTE)
2064 {
2065 vil_image_view<unsigned char > temp = image.get_view();
2066 return temp;
2067 }
2068
2069 if (image.nplanes()==1&&image.pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
2070 {
2071 vil_image_view<unsigned short > temp = image.get_view();
2072 unsigned short vmin=0, vmax=255;
2073 vil_math_value_range<unsigned short>(temp, vmin, vmax);
2074 return brip_vil_float_ops::convert_to_byte(temp, vmin, vmax);
2075 }
2076 // the image is color so we should convert it to greyscale
2077 // Here we assume the color elements are unsigned char.
2078 if (image.nplanes()==3&&image.pixel_format()==VIL_PIXEL_FORMAT_BYTE)
2079 {
2080 vil_image_view<vil_rgb<vxl_byte> > color_image = image.get_view();
2081 unsigned width = color_image.ni(), height = color_image.nj();
2082 // the output image
2083 vil_image_view<unsigned char> grey_image;
2084 grey_image.set_size(width, height);
2085 for (unsigned y = 0; y<height; y++)
2086 for (unsigned x = 0; x<width; x++)
2087 grey_image(x,y) = color_image(x,y).grey();
2088 return grey_image;
2089 }
2090 if (image.nplanes()==3&&image.pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
2091 {
2092 vil_image_view<float> color_image = image.get_view();
2093 unsigned width = color_image.ni(), height = color_image.nj();
2094 // the output image
2095 vil_image_view<float> grey_image_f;
2096 grey_image_f.set_size(width, height);
2097 for (unsigned y = 0; y<height; y++)
2098 for (unsigned x = 0; x<width; x++) {
2099 float v = 0;
2100 for (unsigned p = 0; p<3; ++p)
2101 v += color_image(x,y,p);
2102 grey_image_f(x,y) = v/3.0f;
2103 }
2104 float vmin=0.f, vmax=256.f;
2105 vil_math_value_range<float>(grey_image_f, vmin, vmax);
2106 return brip_vil_float_ops::convert_to_byte(grey_image_f, vmin, vmax);
2107 }
2108 // If we get here then the input is not a type we handle so return a null view
2109 return vil_image_view<vxl_byte>();
2110 }
2111
2112 //--------------------------------------------------------------
2113 //: Read a convolution kernel from file
2114 // Assumes a square kernel with odd dimensions, i.e., w,h = 2n+1
2115 // format:
2116 // \verbatim
2117 // n
2118 // scale
2119 // k00 k01 ... k02n
2120 // ...
2121 // k2n0 k2n1 ... k2n2n
2122 // \endverbatim
2123 //
load_kernel(std::string const & file)2124 vbl_array_2d<float> brip_vil_float_ops::load_kernel(std::string const& file)
2125 {
2126 std::ifstream instr(file.c_str(), std::ios::in);
2127 if (!instr)
2128 {
2129 std::cout << "In brip_vil_float_ops::load_kernel - failed to load kernel\n";
2130 return vbl_array_2d<float>(0,0);
2131 }
2132 unsigned n;
2133 float scale;
2134 float v =0.0f;
2135 instr >> n;
2136 instr >> scale;
2137 unsigned N = 2*n+1;
2138 vbl_array_2d<float> output(N, N);
2139 for (unsigned y = 0; y<N; y++)
2140 for (unsigned x = 0; x<N; x++)
2141 {
2142 instr >> v;
2143 output.put(x, y, v/scale);
2144 }
2145 #ifdef DEBUG
2146 std::cout << "The Kernel\n";
2147 for (unsigned y = 0; y<N; y++)
2148 {
2149 for (unsigned x = 0; x<N; x++)
2150 std::cout << ' ' << output[x][y];
2151 std::cout << '\n';
2152 }
2153 #endif
2154 return output;
2155 }
2156
insert_image(vil_image_view<float> const & image,int col,vnl_matrix<float> & I)2157 static void insert_image(vil_image_view<float> const& image, int col,
2158 vnl_matrix<float> & I)
2159 {
2160 unsigned width = image.ni(), height = image.nj(), row=0;
2161 for (unsigned y =0; y<height; y++)
2162 for (unsigned x = 0; x<width; x++, row++)
2163 I.put(row, col, image(x,y));
2164 }
2165
2166 void brip_vil_float_ops::
basis_images(std::vector<vil_image_view<float>> const & input_images,std::vector<vil_image_view<float>> & basis)2167 basis_images(std::vector<vil_image_view<float> > const& input_images,
2168 std::vector<vil_image_view<float> >& basis)
2169 {
2170 basis.clear();
2171 unsigned n_images = input_images.size();
2172 if (!n_images)
2173 {
2174 std::cout << "In brip_vil_float_ops::basis_images(.) - no input images\n";
2175 return;
2176 }
2177 unsigned width = input_images[0].ni(), height = input_images[0].nj();
2178 unsigned npix = width*height;
2179
2180 // Insert the images into matrix I
2181 vnl_matrix<float> I(npix, n_images, 0.f);
2182 for (unsigned i = 0; i<n_images; i++)
2183 insert_image(input_images[i], i, I);
2184
2185 // Compute the SVD of matrix I
2186 #ifdef DEBUG
2187 std::cout << "Computing Singular values of a " << npix << " by "
2188 << n_images << " matrix\n";
2189 vul_timer t;
2190 #endif
2191 vnl_svd<float> svd(I);
2192 #ifdef DEBUG
2193 std::cout << "SVD Took " << t.real() << " msecs\n"
2194 << "Eigenvalues:\n";
2195 for (unsigned i = 0; i<n_images; i++)
2196 std::cout << svd.W(i) << '\n';
2197 #endif
2198 // Extract the Basis images
2199 unsigned rank = svd.rank();
2200 if (!rank)
2201 {
2202 std::cout << "In brip_vil_float_ops::basis_images(.) - I has zero rank\n";
2203 return;
2204 }
2205 vnl_matrix<float> U = svd.U();
2206 // Output the basis images
2207 unsigned rows = U.rows();
2208 for (unsigned k = 0; k<rank; k++)
2209 {
2210 vil_image_view<float> out(width, height);
2211 unsigned x =0, y = 0;
2212 for (unsigned r = 0; r<rows; r++)
2213 {
2214 out(x, y) = U(r,k);
2215 x++;
2216 if (x>=width)
2217 {
2218 y++;
2219 x=0;
2220 }
2221 if (y>=width)
2222 {
2223 std::cout<<"In brip_vil_float_ops::basis_images(.) - shouldn't happen\n";
2224 return;
2225 }
2226 }
2227 basis.push_back(out);
2228 }
2229 }
2230
2231 //: 1d fourier transform
2232 //-------------------------------------------------------------------------
2233 // This computes an in-place complex-to-complex FFT
2234 // x and y are the real and imaginary arrays of 2^m points.
2235 // dir = 1 gives forward transform
2236 // dir = -1 gives reverse transform
2237 //
2238 // Formula: forward
2239 // \verbatim
2240 // N-1
2241 // ---
2242 // 1 \ - j k 2 pi n / N
2243 // X(n) = --- > x(k) e = forward transform
2244 // N / n=0..N-1
2245 // ---
2246 // k=0
2247 // \endverbatim
2248 //
2249 // Formula: reverse
2250 // \verbatim
2251 // N-1
2252 // ---
2253 // \ j k 2 pi n / N
2254 // X(n) = > x(k) e = forward transform
2255 // / n=0..N-1
2256 // ---
2257 // k=0
2258 // \endverbatim
2259 //
fft_1d(int dir,int m,double * x,double * y)2260 bool brip_vil_float_ops::fft_1d(int dir, int m, double* x, double* y)
2261 {
2262 long nn,i,i1,j,k,i2,l,l1,l2;
2263 double c1,c2,tx,ty,t1,t2,u1,u2,z;
2264
2265 // Calculate the number of points
2266 nn = 1;
2267 for (i=0;i<m;i++)
2268 nn *= 2;
2269
2270 // Do the bit reversal
2271 i2 = nn >> 1;
2272 j = 0;
2273 for (i=0;i<nn-1;i++) {
2274 if (i < j) {
2275 tx = x[i];
2276 ty = y[i];
2277 x[i] = x[j];
2278 y[i] = y[j];
2279 x[j] = tx;
2280 y[j] = ty;
2281 }
2282 k = i2;
2283 while (k <= j) {
2284 j -= k;
2285 k >>= 1;
2286 }
2287 j += k;
2288 }
2289
2290 // Compute the FFT
2291 c1 = -1.0;
2292 c2 = 0.0;
2293 l2 = 1;
2294 for (l=0;l<m;l++) {
2295 l1 = l2;
2296 l2 <<= 1;
2297 u1 = 1.0;
2298 u2 = 0.0;
2299 for (j=0;j<l1;j++) {
2300 for (i=j;i<nn;i+=l2) {
2301 i1 = i + l1;
2302 t1 = u1 * x[i1] - u2 * y[i1];
2303 t2 = u1 * y[i1] + u2 * x[i1];
2304 x[i1] = x[i] - t1;
2305 y[i1] = y[i] - t2;
2306 x[i] += t1;
2307 y[i] += t2;
2308 }
2309 z = u1 * c1 - u2 * c2;
2310 u2 = u1 * c2 + u2 * c1;
2311 u1 = z;
2312 }
2313 c2 = std::sqrt((1.0 - c1) / 2.0);
2314 if (dir == 1)
2315 c2 = -c2;
2316 c1 = std::sqrt((1.0 + c1) / 2.0);
2317 }
2318
2319 // Scaling for forward transform
2320 if (dir == 1) {
2321 for (i=0;i<nn;i++) {
2322 x[i] /= (double)nn;
2323 y[i] /= (double)nn;
2324 }
2325 }
2326
2327 return true;
2328 }
2329
2330 //-------------------------------------------------------------------------
2331 //: Perform a 2D FFT inplace given a complex 2D array
2332 // The direction dir, 1 for forward, -1 for reverse
2333 // The size of the array (nx,ny)
2334 // Return false if there are memory problems or
2335 // the dimensions are not powers of 2
2336 //
fft_2d(vnl_matrix<std::complex<double>> & c,int nx,int ny,int dir)2337 bool brip_vil_float_ops::fft_2d(vnl_matrix<std::complex<double> >& c,int nx,int ny,int dir)
2338 {
2339 int i,j;
2340 int mx, my;
2341 double *real,*imag;
2342 vnl_fft_prime_factors<double> pfx (nx);
2343 vnl_fft_prime_factors<double> pfy (ny);
2344 mx = (int)pfx.pqr()[0];
2345 my = (int)pfy.pqr()[0];
2346 // Transform the rows
2347 real = new double[nx];
2348 imag = new double[nx];
2349 if (real == nullptr || imag == nullptr)
2350 return false;
2351 for (j=0;j<ny;j++) {
2352 for (i=0;i<nx;i++) {
2353 real[i] = c[j][i].real();
2354 imag[i] = c[j][i].imag();
2355 }
2356 brip_vil_float_ops::fft_1d(dir,mx,real,imag);
2357 for (i=0;i<nx;i++) {
2358 std::complex<double> v(real[i], imag[i]);
2359 c[j][i] = v;
2360 }
2361 }
2362 delete [] real;
2363 delete [] imag;
2364 // Transform the columns
2365 real = new double[ny];
2366 imag = new double[ny];
2367 if (real == nullptr || imag == nullptr)
2368 return false;
2369 for (i=0;i<nx;i++) {
2370 for (j=0;j<ny;j++) {
2371 real[j] = c[j][i].real();
2372 imag[j] = c[j][i].imag();
2373 }
2374 fft_1d(dir,my,real,imag);
2375 for (j=0;j<ny;j++) {
2376 std::complex<double> v(real[j], imag[j]);
2377 c[j][i] = v;
2378 }
2379 }
2380 delete [] real;
2381 delete [] imag;
2382 return true;
2383 }
2384
2385 //: reorder the transform values to sequential frequencies as in conventional Fourier transforms.
2386 // The transformation is its self-inverse.
2387 void brip_vil_float_ops::
ftt_fourier_2d_reorder(vnl_matrix<std::complex<double>> const & F1,vnl_matrix<std::complex<double>> & F2)2388 ftt_fourier_2d_reorder(vnl_matrix<std::complex<double> > const& F1,
2389 vnl_matrix<std::complex<double> > & F2)
2390 {
2391 int rows = static_cast<int>(F1.rows()), cols = static_cast<int>(F1.cols());
2392 int half_rows = rows/2, half_cols = cols/2;
2393 int ri, ci;
2394 for (int r = 0; r<rows; r++)
2395 {
2396 if (r<half_rows)
2397 ri = half_rows+r;
2398 else
2399 ri = r-half_rows;
2400 for (int c = 0; c<cols; c++)
2401 {
2402 if (c<half_cols)
2403 ci = half_cols+c;
2404 else
2405 ci = c-half_cols;
2406 F2[ri][ci]=F1[r][c];
2407 }
2408 }
2409 }
2410
2411 //: Compute the fourier transform.
2412 // If the image dimensions are not a power of 2 then the operation fails.
2413 bool brip_vil_float_ops::
fourier_transform(vil_image_view<float> const & input,vil_image_view<float> & mag,vil_image_view<float> & phase)2414 fourier_transform(vil_image_view<float> const& input,
2415 vil_image_view<float>& mag,
2416 vil_image_view<float>& phase)
2417 {
2418 unsigned w = input.ni(), h = input.nj();
2419
2420 vnl_fft_prime_factors<float> pfx (w);
2421 vnl_fft_prime_factors<float> pfy (h);
2422 if (!pfx.pqr()[0]||!pfy.pqr()[0])
2423 return false;
2424 // fill the fft matrix
2425 vnl_matrix<std::complex<double> > fft_matrix(h, w), fourier_matrix(h,w);
2426 for (unsigned y = 0; y<h; y++)
2427 for (unsigned x =0; x<w; x++)
2428 {
2429 std::complex<double> cv(input(x,y), 0.0);
2430 fft_matrix.put(y, x, cv);
2431 }
2432 #ifdef DEBUG
2433 for (unsigned r = 0; r<h; r++)
2434 for (unsigned c =0; c<w; c++)
2435 {
2436 std::complex<double> res = fft_matrix[r][c];
2437 std::cout << res << '\n';
2438 }
2439 #endif
2440
2441 brip_vil_float_ops::fft_2d(fft_matrix, w, h, 1);
2442 brip_vil_float_ops::ftt_fourier_2d_reorder(fft_matrix, fourier_matrix);
2443 mag.set_size(w,h);
2444 phase.set_size(w,h);
2445
2446 // extract magnitude and phase
2447 for (unsigned r = 0; r<h; r++)
2448 for (unsigned c = 0; c<w; c++)
2449 {
2450 auto re = (float)fourier_matrix[r][c].real(),
2451 im = (float)fourier_matrix[r][c].imag();
2452 mag(c,r) = std::sqrt(re*re + im*im);
2453 phase(c,r) = std::atan2(im, re);
2454 }
2455 return true;
2456 }
2457
2458 bool brip_vil_float_ops::
inverse_fourier_transform(vil_image_view<float> const & mag,vil_image_view<float> const & phase,vil_image_view<float> & output)2459 inverse_fourier_transform(vil_image_view<float> const& mag,
2460 vil_image_view<float> const& phase,
2461 vil_image_view<float>& output)
2462 {
2463 unsigned w = mag.ni(), h = mag.nj();
2464 vnl_matrix<std::complex<double> > fft_matrix(h, w), fourier_matrix(h, w);
2465 for (unsigned y = 0; y<h; y++)
2466 for (unsigned x =0; x<w; x++)
2467 {
2468 float m = mag(x,y);
2469 float p = phase(x,y);
2470 std::complex<double> cv(m*std::cos(p), m*std::sin(p));
2471 fourier_matrix.put(y, x, cv);
2472 }
2473
2474 brip_vil_float_ops::ftt_fourier_2d_reorder(fourier_matrix, fft_matrix);
2475 brip_vil_float_ops::fft_2d(fft_matrix, w, h, -1);
2476
2477 output.set_size(w,h);
2478
2479 for (unsigned y = 0; y<h; y++)
2480 for (unsigned x = 0; x<w; x++)
2481 output(x,y) = (float)fft_matrix[y][x].real();
2482 return true;
2483 }
2484
resize(vil_image_view<float> const & input,unsigned width,unsigned height,vil_image_view<float> & output)2485 void brip_vil_float_ops::resize(vil_image_view<float> const& input,
2486 unsigned width, unsigned height,
2487 vil_image_view<float>& output)
2488 {
2489 unsigned w = input.ni(), h = input.nj();
2490 output.set_size(width, height);
2491 for (unsigned y = 0; y<height; y++)
2492 for (unsigned x = 0; x<width; x++)
2493 if (x<w && y<h)
2494 output(x,y) = input(x,y);
2495 else
2496 output(x,y) = 0; // pad with zeroes
2497 }
2498
2499 //: resize the input to the closest power of two image dimensions
2500 bool brip_vil_float_ops::
resize_to_power_of_two(vil_image_view<float> const & input,vil_image_view<float> & output)2501 resize_to_power_of_two(vil_image_view<float> const& input,
2502 vil_image_view<float>& output)
2503 {
2504 unsigned max_exp = 13; // we wouldn't want to have such large images in memory
2505 unsigned w = input.ni(), h = input.nj();
2506 unsigned prodw = 1, prodh = 1;
2507 // Find power of two width
2508 unsigned nw, nh;
2509 for (nw = 1; nw<=max_exp; nw++)
2510 if (prodw>w)
2511 break;
2512 else
2513 prodw *= 2;
2514 if (nw==max_exp)
2515 return false;
2516 // Find power of two height
2517 for (nh = 1; nh<=max_exp; nh++)
2518 if (prodh>h)
2519 break;
2520 else
2521 prodh *= 2;
2522 if (nh==max_exp)
2523 return false;
2524 brip_vil_float_ops::resize(input, prodw, prodh, output);
2525
2526 return true;
2527 }
2528
2529 //
2530 //: block a periodic signal by suppressing two Gaussian lobes in the frequency domain.
2531 // The lobes are on the line defined by dir_fx and dir_fy through the
2532 // dc origin, assumed (0, 0). The center frequency, f0, is the distance along
2533 // the line to the center of each blocking lobe (+- f0). radius is the
2534 // standard deviation of each lobe. Later we can define a "filter" class.
2535 //
gaussian_blocking_filter(float dir_fx,float dir_fy,float f0,float radius,float fx,float fy)2536 float brip_vil_float_ops::gaussian_blocking_filter(float dir_fx,
2537 float dir_fy,
2538 float f0,
2539 float radius,
2540 float fx,
2541 float fy)
2542 {
2543 // normalize dir_fx and dir_fy
2544 float mag = std::sqrt(dir_fx*dir_fx + dir_fy*dir_fy);
2545 if (!mag)
2546 return 0.f;
2547 float r2 = 2.f*radius*radius;
2548 float dx = dir_fx/mag, dy = dir_fy/mag;
2549 // compute the centers of each lobe
2550 float fx0p = dx*f0, fy0p = dy*f0;
2551 float fx0m = -dx*f0, fy0m = -dy*f0;
2552 // the squared distance of fx, fy from each center
2553 float d2p = (fx-fx0p)*(fx-fx0p) + (fy-fy0p)*(fy-fy0p);
2554 float d2m = (fx-fx0m)*(fx-fx0m) + (fy-fy0m)*(fy-fy0m);
2555 // use the closest lobe
2556 float d = d2p;
2557 if (d2m<d2p)
2558 d = d2m;
2559 // the gaussian blocking function
2560 float gb = 1.f-(float)std::exp(-d/r2);
2561 return gb;
2562 }
2563
2564 bool brip_vil_float_ops::
spatial_frequency_filter(vil_image_view<float> const & input,float dir_fx,float dir_fy,float f0,float radius,bool output_fourier_mag,vil_image_view<float> & output)2565 spatial_frequency_filter(vil_image_view<float> const& input,
2566 float dir_fx, float dir_fy,
2567 float f0, float radius,
2568 bool output_fourier_mag,
2569 vil_image_view<float>& output)
2570 {
2571 // Compute the fourier transform of the image.
2572 vil_image_view<float> pow_two, mag, bmag, phase, pow_two_filt;
2573 brip_vil_float_ops::resize_to_power_of_two(input, pow_two);
2574 unsigned Nfx = pow_two.ni(), Nfy = pow_two.nj();
2575
2576 if (!brip_vil_float_ops::fourier_transform(pow_two, mag, phase))
2577 return false;
2578 bmag.set_size(Nfx, Nfy);
2579
2580 // filter the magnitude function
2581 float Ofx = Nfx*0.5f, Ofy = Nfy*0.5f;
2582 for (unsigned fy =0; fy<Nfy; fy++)
2583 for (unsigned fx =0; fx<Nfx; fx++)
2584 {
2585 float gb = gaussian_blocking_filter(dir_fx, dir_fy, f0,
2586 radius,
2587 (float)fx-Ofx, (float)fy-Ofy);
2588 bmag(fx,fy) = mag(fx,fy)*gb;
2589 }
2590 if (output_fourier_mag)
2591 {
2592 output = bmag;
2593 return true;
2594 }
2595 // Transform back
2596 pow_two_filt.set_size(Nfx, Nfy);
2597 brip_vil_float_ops::inverse_fourier_transform(bmag, phase, pow_two_filt);
2598
2599 // Resize to original input size
2600 brip_vil_float_ops::resize(pow_two_filt, input.ni(), input.nj(), output);
2601 return true;
2602 }
2603
2604 //----------------------------------------------------------------------
2605 //: Bi-linear interpolation on the neighborhood below.
2606 // \verbatim
2607 // xr
2608 // yr 0 x
2609 // x x
2610 // \endverbatim
2611 //
2612 double brip_vil_float_ops::
bilinear_interpolation(vil_image_view<float> const & input,double x,double y)2613 bilinear_interpolation(vil_image_view<float> const& input, double x, double y)
2614 {
2615 // check bounds
2616 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
2617 // the pixel containing the interpolated point
2618 int xr = (int)x, yr = (int)y;
2619 double fx = x-xr, fy = y-yr;
2620 if (xr<0||xr>w-2)
2621 return 0.0;
2622 if (yr<0||yr>h-2)
2623 return 0.0;
2624 double int00 = input(xr, yr), int10 = input(xr+1,yr);
2625 double int01 = input(xr, yr+1), int11 = input(xr+1,yr+1);
2626 double int0 = int00 + fy * (int01 - int00);
2627 double int1 = int10 + fy * (int11 - int10);
2628 auto val = (float) (int0 + fx * (int1 - int0));
2629 return val;
2630 }
2631
2632 //: Transform the input to the output by a homography.
2633 // If the output size is fixed then only the corresponding
2634 // region of input image space is transformed.
homography(vil_image_view<float> const & input,vgl_h_matrix_2d<double> const & H,vil_image_view<float> & output,bool output_size_fixed,float output_fill_value)2635 bool brip_vil_float_ops::homography(vil_image_view<float> const& input,
2636 vgl_h_matrix_2d<double> const& H,
2637 vil_image_view<float>& output,
2638 bool output_size_fixed,
2639 float output_fill_value)
2640 {
2641 if (!input)
2642 return false;
2643 // smooth the input to condition interpolation
2644 vil_image_view<float> gimage =
2645 brip_vil_float_ops::gaussian(input, 0.5f);
2646
2647 // First, there is some rather complex bookeeping to insure that
2648 // the input and output image rois are consistent with the homography.
2649
2650 // the bounding boxes corresponding to input and output rois
2651 // We also construct polygons since homographies turn boxes into arbitrary
2652 // quadrilaterals.
2653 vsol_box_2d_sptr input_roi, output_roi;
2654 vsol_polygon_2d_sptr input_poly, output_poly;
2655 vgl_h_matrix_2d<double> Hinv;
2656 // set up the roi and poly for the input image
2657 unsigned win = gimage.ni(), hin = gimage.nj();
2658 input_roi = new vsol_box_2d();
2659 input_roi->add_point(0, 0);
2660 input_roi->add_point(win, hin);
2661 input_poly = bsol_algs::poly_from_box(input_roi);
2662 // Case I
2663 // the output image size and input transform can be adjusted
2664 // to map the transformed image onto the full range
2665 if (!output_size_fixed)
2666 {
2667 if (!bsol_algs::homography(input_poly, H, output_poly))
2668 return false;
2669 vsol_box_2d_sptr temp = output_poly->get_bounding_box();
2670 output.set_size((int)temp->width(), (int)temp->height());
2671 output.fill(output_fill_value);
2672 //offset the transform so the origin is (0,0)
2673 output_roi = new vsol_box_2d();
2674 output_roi->add_point(0, 0);
2675 output_roi->add_point(temp->width(), temp->height());
2676 vgl_h_matrix_2d<double> t; t.set_identity().set_translation(-temp->get_min_x(),-temp->get_min_y());
2677 Hinv = (t*H).get_inverse();
2678 }
2679 else // Case II, the output image size is fixed so we have to find the
2680 { // inverse mapping of the output roi and intersect with the input roi
2681 // to determine the domain of the mapping
2682 if (!output)
2683 return false;
2684 //The output roi and poly
2685 unsigned wout = output.ni(), hout = output.nj();
2686 output.fill(output_fill_value);
2687 output_roi = new vsol_box_2d();
2688 output_roi->add_point(0, 0);
2689 output_roi->add_point(wout, hout);
2690 output_poly = bsol_algs::poly_from_box(output_roi);
2691
2692 //Construct the reverse mapping of the output bounds
2693 vsol_polygon_2d_sptr tpoly;
2694 Hinv = H.get_inverse();
2695 if (!bsol_algs::homography(output_poly, Hinv, tpoly))
2696 return false;
2697
2698 //form the roi corresponding to the inverse mapped output bounds
2699 vsol_box_2d_sptr tbox = tpoly->get_bounding_box();
2700
2701 //intersect with the input image bounds to get the input roi
2702 vsol_box_2d_sptr temp;
2703 if (!bsol_algs::intersection(tbox, input_roi, temp))
2704 return false;
2705 input_roi = temp;
2706 }
2707 // At this point we have the correct bounds for the input and
2708 // the output image
2709
2710 // Iterate over the output image space and map the location of each
2711 // pixel into the input image space. Then carry out interpolation to
2712 // get the value of each output pixel
2713
2714 // Dimensions of the input image
2715 int ailow = int(input_roi->get_min_x()+0.9999f); // round up to nearest int
2716 int aihigh = int(input_roi->get_max_x()); // round down to nearest int
2717 int ajlow = int(input_roi->get_min_y()+0.9999f);
2718 int ajhigh = int(input_roi->get_max_y());
2719
2720 // Dimensions of the output image
2721 int bilow = int(output_roi->get_min_x()+0.9999f);
2722 int bihigh = int(output_roi->get_max_x());
2723 int bjlow = int(output_roi->get_min_y()+0.9999f);
2724 int bjhigh = int(output_roi->get_max_y());
2725
2726 // The inverse transform is used to map backwards from the output
2727 const vnl_matrix_fixed<double,3,3>& Minv = Hinv.get_matrix();
2728
2729 // Now use Hinv to transform the image
2730 for (int i = bilow; i<bihigh; i++)
2731 for (int j = bjlow; j<bjhigh; j++)
2732 {
2733 // Transform the pixel
2734 float val;
2735 double u = Minv[0][0] * i + Minv[0][1] * j + Minv[0][2];
2736 double v = Minv[1][0] * i + Minv[1][1] * j + Minv[1][2];
2737 double w = Minv[2][0] * i + Minv[2][1] * j + Minv[2][2];
2738 u /= w;
2739 v /= w;
2740
2741 // Now do linear interpolation
2742 {
2743 int iu = (int) u;
2744 int iv = (int) v;
2745
2746 if (iu >= ailow && iu < aihigh-1 &&
2747 iv >= ajlow && iv < ajhigh-1)
2748 {
2749 // Get the neighbouring pixels
2750 // (u v) (u+1 v)
2751 // (u v+1) (u+1 v+1)
2752 //
2753 double v00 = gimage(iu, iv);
2754 double v01 = gimage(iu, iv+1);
2755 double v10 = gimage(iu+1,iv);
2756 double v11 = gimage(iu+1, iv+1);
2757
2758 double fu = u - iu;
2759 double fv = v - iv;
2760 double v0 = v00 + fv * (v01 - v00);
2761 double v1 = v10 + fv * (v11 - v10);
2762 val = (float) (v0 + fu * (v1 - v0));
2763 // Set the value
2764 output(i,j) = val;
2765 }
2766 }
2767 }
2768 return true;
2769 }
2770
2771 //: rotate the input image counter-clockwise about the image origin.
2772 // Demonstrates the use of image homography
2773 vil_image_view<float>
rotate(vil_image_view<float> const & input,double theta_deg)2774 brip_vil_float_ops::rotate(vil_image_view<float> const& input,
2775 double theta_deg)
2776 {
2777 vil_image_view<float> out;
2778 if (!input)
2779 return out;
2780 double ang = theta_deg;
2781 // map theta_deg to [0 360]
2782 while (ang>360)
2783 ang-=360;
2784 while (ang<0)
2785 ang+=360;
2786 // convert to radians
2787 double deg_to_rad = vnl_math::pi_over_180;
2788 double rang = deg_to_rad*ang;
2789 vgl_h_matrix_2d<double> H;
2790 H.set_identity().set_rotation(rang);
2791 vil_image_view<float> temp;
2792 // The transform is adjusted to map the full input domain onto
2793 // the output image.
2794 if (!brip_vil_float_ops::homography(input, H, temp))
2795 return out;
2796 return temp;
2797 }
2798
chip(vil_image_view<float> const & input,vsol_box_2d_sptr const & roi,vil_image_view<float> & chip)2799 bool brip_vil_float_ops::chip(vil_image_view<float> const& input,
2800 vsol_box_2d_sptr const& roi,
2801 vil_image_view<float>& chip)
2802 {
2803 if (!input||!roi)
2804 return false;
2805 int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
2806 int x_min = (int)roi->get_min_x(), y_min = (int)roi->get_min_y();
2807 int x_max = (int)roi->get_max_x(), y_max = (int)roi->get_max_y();
2808 if (x_min<0)
2809 x_min = 0;
2810 if (y_min<0)
2811 y_min = 0;
2812 if (x_max>w-1)
2813 x_max=w-1;
2814 if (y_max>h-1)
2815 y_max=h-1;
2816 int rw = x_max-x_min, rh = y_max-y_min;
2817 if (rw<=0||rh<=0)
2818 return false;
2819 vil_image_view<float> temp(rw+1, rh+1, 1);
2820 for (int y = y_min; y<=y_max; y++) //!< changed < to <= to include the boundary points too
2821 for (int x =x_min; x<=x_max; x++)
2822 temp(x-x_min, y-y_min) = input(x, y);
2823 chip = temp;
2824 return true;
2825 }
2826
2827 //: convert image resource to cropped view according to a roi.
chip(vil_image_resource_sptr const & image,brip_roi_sptr const & roi,vil_image_resource_sptr & chip)2828 bool brip_vil_float_ops::chip(vil_image_resource_sptr const& image,
2829 brip_roi_sptr const& roi,
2830 vil_image_resource_sptr & chip)
2831 {
2832 // image bounds
2833 unsigned ni = image->ni(), nj = image->nj();
2834
2835 // image bounds for the chip
2836 unsigned cm = roi->cmin(0), rm = roi->rmin(0);
2837 unsigned niv = roi->csize(0), njv = roi->rsize(0);
2838 std::cout << "org(" << cm << ' ' << rm << "):size(" << niv
2839 << ' ' << njv << ")\n" << std::flush;
2840 // check bounds
2841 if (cm>ni-1||rm>nj-1||cm+niv>ni||rm+njv>nj)
2842 return false;
2843 vil_pixel_format pix_format = image->pixel_format();
2844 // get an appropriate image view for scalar images we care about
2845 if (image->nplanes()==1)
2846 {
2847 if (pix_format==VIL_PIXEL_FORMAT_BYTE)
2848 {
2849 vil_image_view<unsigned char> temp = image->get_view(cm, niv, rm, njv);
2850 if (!temp) return false;
2851 chip = vil_new_image_resource_of_view(temp);
2852 return true;
2853 }
2854 else if (pix_format==VIL_PIXEL_FORMAT_UINT_16)
2855 {
2856 vil_image_view<unsigned short> temp = image->get_view(cm, niv, rm, njv);
2857 if (!temp) return false;
2858 chip = vil_new_image_resource_of_view(temp);
2859 std::cout << "resc(" << chip->ni() << ' ' << chip->nj()<< ")\n"
2860 << std::flush;
2861 return true;
2862 }
2863 else if (pix_format==VIL_PIXEL_FORMAT_FLOAT)
2864 {
2865 vil_image_view<float> temp = image->get_view(cm, niv, rm, njv);
2866 if (!temp) return false;
2867 chip = vil_new_image_resource_of_view(temp);
2868 return true;
2869 }
2870 }
2871
2872 // color data
2873 if (image->nplanes()==3)
2874 {
2875 if (pix_format==VIL_PIXEL_FORMAT_BYTE) //the only way now
2876 {
2877 //extract view corresponding to region of interest
2878 vil_image_view<vil_rgb<vxl_byte> > temp = image->get_view(cm, niv, rm, njv);
2879 if (!temp) return false;
2880 chip = vil_new_image_resource_of_view(temp);
2881 return true;
2882 }
2883 }
2884 return false;
2885 }
2886
2887 bool brip_vil_float_ops::
chip(std::vector<vil_image_resource_sptr> const & images,brip_roi_sptr const & roi,std::vector<vil_image_resource_sptr> & chips)2888 chip(std::vector<vil_image_resource_sptr> const& images,
2889 brip_roi_sptr const& roi,
2890 std::vector<vil_image_resource_sptr>& chips)
2891 {
2892 for (const auto & image : images) {
2893 vil_image_resource_sptr chip;
2894 if (!brip_vil_float_ops::chip(image, roi, chip))
2895 return false;
2896 chips.push_back(chip);
2897 }
2898 return true;
2899 }
2900
2901 //: perform normalized cross-correlation at a sub-pixel location.
2902 // Thus all the pixel values are interpolated.
2903 float brip_vil_float_ops::
cross_correlate(vil_image_view<float> const & image1,vil_image_view<float> const & image2,float x,float y,int radius,float intensity_thresh)2904 cross_correlate(vil_image_view<float> const& image1,
2905 vil_image_view<float> const& image2,
2906 float x, float y, int radius, float intensity_thresh)
2907 {
2908 unsigned w1 = image1.ni(), h1 = image1.nj();
2909 unsigned w2 = image1.ni(), h2 = image1.nj();
2910 // bounds checks
2911 if (w1!=w2||h1!=h2)
2912 return -1;
2913 if (x<radius||x>w1-radius-1||y<radius||y>h1-radius-1)
2914 return -1;
2915
2916 // accumulate correlation sums,
2917 // bi-linear interpolate the values
2918 int s = 2*radius+1;
2919 double area = s*s;
2920 double sI1=0, sI2=0, sI1I1=0, sI2I2=0, sI1I2=0;
2921 for (int y0 = -10*radius; y0<=10*radius; ++y0)
2922 for (int x0 = -10*radius; x0<=10*radius; ++x0)
2923 {
2924 float xp = x+0.1f*(float)x0, yp = y+0.1f*(float)y0;
2925 double v1 = brip_vil_float_ops::bilinear_interpolation(image1, xp, yp);
2926 double v2 = brip_vil_float_ops::bilinear_interpolation(image2, xp, yp);
2927 sI1 += v1;
2928 sI2 += v2;
2929 sI1I1 += v1*v1;
2930 sI2I2 += v2*v2;
2931 sI1I2 += v1*v2;
2932 }
2933 // compute correlation.
2934 float cc = cross_corr(area, sI1, sI2, sI1I1, sI2I2, sI1I2, intensity_thresh);
2935 return cc;
2936 }
2937
2938 //: r0 is the image from which to read the new intensity values.
2939 // \param r is the summing array row in which the values are to be accumulated
update_row(vil_image_view<float> const & image1,vil_image_view<float> const & image2,int r0,int r,vbl_array_2d<double> & SI1,vbl_array_2d<double> & SI2,vbl_array_2d<double> & SI1I1,vbl_array_2d<double> & SI2I2,vbl_array_2d<double> & SI1I2)2940 static bool update_row(vil_image_view<float> const& image1,
2941 vil_image_view<float> const& image2,
2942 int r0, int r,
2943 vbl_array_2d<double>& SI1,
2944 vbl_array_2d<double>& SI2,
2945 vbl_array_2d<double>& SI1I1,
2946 vbl_array_2d<double>& SI2I2,
2947 vbl_array_2d<double>& SI1I2)
2948 {
2949 unsigned w1 = image1.ni();
2950 unsigned w2 = image2.ni();
2951 int h1 = image1.nj();
2952 int h2 = image2.nj();
2953 if (w1!=w2||h1!=h2||r<0||r>=h1)
2954 return false;
2955 double i10 = image1(0,r0), i20 = image2(0,r0);
2956 SI1[r][0] = i10; SI2[r][0] = i20; SI1I1[r][0]=i10*i10;
2957 SI2I2[r][0]=i20*i20; SI1I2[r][0]=i10*i20;
2958 for (unsigned c = 1; c<w1; c++)
2959 {
2960 double i1c = image1(c,r0);
2961 double i2c = image2(c,r0);
2962 SI1[r][c] = SI1[r][c-1]+i1c;
2963 SI2[r][c] = SI2[r][c-1]+i2c;
2964 SI1I1[r][c] = SI1I1[r][c-1]+ i1c*i1c;
2965 SI2I2[r][c] = SI2I2[r][c-1]+ i2c*i2c;
2966 SI1I2[r][c] = SI1I2[r][c-1]+ i1c*i2c;
2967 }
2968 return true;
2969 }
2970
initialize_slice(vil_image_view<float> const & image1,vil_image_view<float> const & image2,unsigned radius,vbl_array_2d<double> & SI1,vbl_array_2d<double> & SI2,vbl_array_2d<double> & SI1I1,vbl_array_2d<double> & SI2I2,vbl_array_2d<double> & SI1I2)2971 static bool initialize_slice(vil_image_view<float> const& image1,
2972 vil_image_view<float> const& image2,
2973 unsigned radius,
2974 vbl_array_2d<double>& SI1,
2975 vbl_array_2d<double>& SI2,
2976 vbl_array_2d<double>& SI1I1,
2977 vbl_array_2d<double>& SI2I2,
2978 vbl_array_2d<double>& SI1I2)
2979 {
2980 for (unsigned r = 0; r<=2*radius; r++)
2981 if (!update_row(image1, image2, r, r, SI1, SI2, SI1I1, SI2I2, SI1I2))
2982 return false;
2983 return true;
2984 }
2985
collapse_slice(vbl_array_2d<double> const & SI1,vbl_array_2d<double> const & SI2,vbl_array_2d<double> const & SI1I1,vbl_array_2d<double> const & SI2I2,vbl_array_2d<double> const & SI1I2,vbl_array_1d<double> & dSI1,vbl_array_1d<double> & dSI2,vbl_array_1d<double> & dSI1I1,vbl_array_1d<double> & dSI2I2,vbl_array_1d<double> & dSI1I2)2986 static bool collapse_slice(vbl_array_2d<double> const& SI1,
2987 vbl_array_2d<double> const& SI2,
2988 vbl_array_2d<double> const& SI1I1,
2989 vbl_array_2d<double> const& SI2I2,
2990 vbl_array_2d<double> const& SI1I2,
2991 vbl_array_1d<double>& dSI1,
2992 vbl_array_1d<double>& dSI2,
2993 vbl_array_1d<double>& dSI1I1,
2994 vbl_array_1d<double>& dSI2I2,
2995 vbl_array_1d<double>& dSI1I2)
2996 {
2997 // sanity check
2998 unsigned w = SI1.cols(), h = SI1.rows();
2999 unsigned dw = dSI1.size();
3000 if (dw!=w)
3001 return false;
3002
3003 for (unsigned c = 0; c<w; c++)
3004 {
3005 dSI1[c]=0; dSI2[c]=0; dSI1I1[c]=0;
3006 dSI2I2[c]=0; dSI1I2[c]=0;
3007 for (unsigned r = 0; r<h; r++)
3008 {
3009 dSI1[c] += SI1[r][c];
3010 dSI2[c] += SI2[r][c];
3011 dSI1I1[c] += SI1I1[r][c];
3012 dSI2I2[c] += SI2I2[r][c];
3013 dSI1I2[c] += SI1I2[r][c];
3014 }
3015 }
3016 return true;
3017 }
3018
cross_correlate_row(int radius,vbl_array_1d<double> & dSI1,vbl_array_1d<double> & dSI2,vbl_array_1d<double> & dSI1I1,vbl_array_1d<double> & dSI2I2,vbl_array_1d<double> & dSI1I2,float intensity_thresh,vbl_array_1d<float> & cc)3019 static bool cross_correlate_row(int radius,
3020 vbl_array_1d<double>& dSI1,
3021 vbl_array_1d<double>& dSI2,
3022 vbl_array_1d<double>& dSI1I1,
3023 vbl_array_1d<double>& dSI2I2,
3024 vbl_array_1d<double>& dSI1I2,
3025 float intensity_thresh,
3026 vbl_array_1d<float>& cc)
3027 {
3028 // sanity check
3029 unsigned w = dSI1.size(), wc = cc.size();
3030 if (!w||!wc||w!=wc)
3031 return false;
3032 int s = 2*radius+1;
3033 double area = s*s;
3034 // the general case
3035 double si1=dSI1[s-1], si2=dSI2[s-1], si1i1=dSI1I1[s-1],
3036 si2i2=dSI2I2[s-1], si1i2=dSI1I2[s-1];
3037 cc[radius]= cross_corr(area, si1, si2, si1i1, si2i2, si1i2, intensity_thresh);
3038 // the remaining columns
3039 for (unsigned c = radius+1; c+radius<w; c++)
3040 {
3041 si1=dSI1[c+radius]-dSI1[c-radius-1];
3042 si2=dSI2[c+radius]-dSI2[c-radius-1];
3043 si1i1=dSI1I1[c+radius]-dSI1I1[c-radius-1];
3044 si2i2=dSI2I2[c+radius]-dSI2I2[c-radius-1];
3045 si1i2=dSI1I2[c+radius]-dSI1I2[c-radius-1];
3046 cc[c] = cross_corr(area, si1, si2, si1i1, si2i2, si1i2, intensity_thresh);
3047 }
3048 return true;
3049 }
3050
advance_rows(vbl_array_2d<double> & S)3051 static void advance_rows(vbl_array_2d<double>& S)
3052 {
3053 unsigned nr = S.rows(), nc = S.cols();
3054 for (unsigned r = 0; r<nr-1; r++)
3055 for (unsigned c =0; c<nc; c++)
3056 S[r][c]=S[r+1][c];
3057 }
3058
output_cc_row(int r0,vbl_array_1d<float> const & cc,vil_image_view<float> & out)3059 static bool output_cc_row(int r0, vbl_array_1d<float> const& cc,
3060 vil_image_view<float>& out)
3061 {
3062 unsigned n = cc.size(), w = out.ni();
3063 if (n!=w)
3064 return false;
3065 for (unsigned c = 0; c<w; c++)
3066 out(c, r0) = cc[c];
3067 return true;
3068 }
3069
3070 bool brip_vil_float_ops::
cross_correlate(vil_image_view<float> const & image1,vil_image_view<float> const & image2,vil_image_view<float> & out,int radius,float intensity_thresh)3071 cross_correlate(vil_image_view<float> const& image1,
3072 vil_image_view<float> const& image2,
3073 vil_image_view<float>& out,
3074 int radius, float intensity_thresh)
3075 {
3076 vul_timer t;
3077 unsigned w = image1.ni(), h = image1.nj();
3078 unsigned w2 = image2.ni(), h2 = image2.nj();
3079 // sizes must match
3080 if (w!=w2||h!=h2)
3081 {
3082 std::cout << "In brip_vil_float_ops::cross_correlate(..) -"
3083 << " image sizes don't match\n";
3084 return out;
3085 }
3086 out.set_size(w, h);
3087 out.fill(0.f);
3088 int s = 2*radius+1;
3089 // Create the running sum slices
3090 vbl_array_2d<double> SI1(s,w), SI2(s,w),
3091 SI1I1(s,w), SI2I2(s,w), SI1I2(s,w);
3092 vbl_array_1d<float> cc(w, 0.f);
3093 vbl_array_1d<double> dSI1(w, 0.0), dSI2(w, 0.0),
3094 dSI1I1(w, 0.0), dSI2I2(w, 0.0), dSI1I2(w, 0.0);
3095 initialize_slice(image1, image2, radius, SI1, SI2, SI1I1, SI2I2, SI1I2);
3096 if (!collapse_slice(SI1, SI2, SI1I1, SI2I2, SI1I2,
3097 dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2))
3098 return false;
3099 unsigned r0 = radius;
3100 for (; r0+radius+1<h; r0++)
3101 {
3102 if (r0==5)
3103 r0=5;
3104 #ifdef DEBUG
3105 std::cout << "r0 " << r0 << '\n';
3106 #endif
3107 if (!cross_correlate_row(radius, dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2,
3108 intensity_thresh, cc))
3109 return false;
3110 #ifdef DEBUG
3111 std::cout << '\n';
3112 #endif
3113 advance_rows(SI1); advance_rows(SI2); advance_rows(SI1I1);
3114 advance_rows(SI2I2); advance_rows(SI1I2);
3115 if (!update_row(image1, image2, r0+radius+1, 2*radius,
3116 SI1, SI2, SI1I1, SI2I2, SI1I2))
3117 return false;
3118 if (!collapse_slice(SI1, SI2, SI1I1, SI2I2, SI1I2,
3119 dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2))
3120 return false;
3121 if (!output_cc_row(r0, cc, out))
3122 return out;
3123 }
3124 // handle the last row
3125 #ifdef DEBUG
3126 std::cout << "r0 " << r0 << '\n';
3127 #endif
3128 if (!cross_correlate_row(radius, dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2,
3129 intensity_thresh, cc))
3130 return false;
3131 #ifdef DEBUG
3132 std::cout << '\n';
3133 #endif
3134 if (!output_cc_row(r0, cc, out))
3135 return false;
3136 #ifdef DEBUG
3137 std::cout << "RunningSumCrossCorrelation for " << w*h/1000.0f << " k pixels in "
3138 << t.real() << " msecs\n"<< std::flush;
3139 #endif
3140 return true;
3141 }
3142
3143 vil_image_resource_sptr
sum(vil_image_resource_sptr const & img0,vil_image_resource_sptr const & img1)3144 brip_vil_float_ops::sum(vil_image_resource_sptr const& img0,
3145 vil_image_resource_sptr const& img1)
3146 {
3147 vil_image_view<float> op0 = brip_vil_float_ops::convert_to_float(img0);
3148 vil_image_view<float> op1 = brip_vil_float_ops::convert_to_float(img1);
3149 vil_image_view<float> sum;
3150 vil_math_image_sum(op0, op1, sum);
3151
3152 // find out the types of the input images for now, only do greyscale operands
3153 vil_pixel_format pix_format0 = img0->pixel_format();
3154 vil_pixel_format pix_format1 = img1->pixel_format();
3155 if (pix_format0 == VIL_PIXEL_FORMAT_FLOAT ||
3156 pix_format1 == VIL_PIXEL_FORMAT_FLOAT)
3157 {
3158 return vil_new_image_resource_of_view(sum);
3159 }
3160
3161 if (pix_format0 == VIL_PIXEL_FORMAT_UINT_16 ||
3162 pix_format1 == VIL_PIXEL_FORMAT_UINT_16)
3163 {
3164 vil_image_view<vxl_uint_16> res =
3165 brip_vil_float_ops::convert_to_short(vil_new_image_resource_of_view(sum));
3166 return vil_new_image_resource_of_view(res);
3167 }
3168
3169 if (pix_format0 == VIL_PIXEL_FORMAT_BYTE ||
3170 pix_format1 == VIL_PIXEL_FORMAT_BYTE)
3171 {
3172 vil_image_view<vxl_byte> res =
3173 brip_vil_float_ops::convert_to_byte(sum);
3174 return vil_new_image_resource_of_view(res);
3175 }
3176 // for color return the float image
3177 return vil_new_image_resource_of_view(sum);
3178 }
3179
3180 // Compute img0 - img1
3181 vil_image_resource_sptr
difference(vil_image_resource_sptr const & img0,vil_image_resource_sptr const & img1)3182 brip_vil_float_ops::difference(vil_image_resource_sptr const& img0,
3183 vil_image_resource_sptr const& img1)
3184 {
3185 vil_image_view<float> op0 = brip_vil_float_ops::convert_to_float(img0);
3186 vil_image_view<float> op1 = brip_vil_float_ops::convert_to_float(img1);
3187 vil_image_view<float> diff;
3188 vil_math_image_difference(op0, op1, diff);
3189
3190 // find out the types of the input images for now, only do greyscale operands
3191 vil_pixel_format pix_format0 = img0->pixel_format();
3192 vil_pixel_format pix_format1 = img1->pixel_format();
3193 if (pix_format0 == VIL_PIXEL_FORMAT_FLOAT ||
3194 pix_format1 == VIL_PIXEL_FORMAT_FLOAT)
3195 {
3196 return vil_new_image_resource_of_view(diff);
3197 }
3198
3199 if (pix_format0 == VIL_PIXEL_FORMAT_UINT_16 ||
3200 pix_format1 == VIL_PIXEL_FORMAT_UINT_16)
3201 {
3202 vil_image_view<vxl_uint_16> res =
3203 brip_vil_float_ops::convert_to_short(vil_new_image_resource_of_view(diff));
3204 return vil_new_image_resource_of_view(res);
3205 }
3206
3207 if (pix_format0 == VIL_PIXEL_FORMAT_BYTE ||
3208 pix_format1 == VIL_PIXEL_FORMAT_BYTE)
3209 {
3210 vil_image_view<vxl_byte> res = brip_vil_float_ops::convert_to_byte(diff);
3211 return vil_new_image_resource_of_view(res);
3212 }
3213 // for color return the float image
3214 return vil_new_image_resource_of_view(diff);
3215 }
3216
3217 //: Compute the entropy of the intensity of a region
3218 // Note no bounds checking!
entropy_i(const unsigned i,const unsigned j,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & intensity,const float range,const unsigned bins)3219 float brip_vil_float_ops::entropy_i(const unsigned i, const unsigned j,
3220 const unsigned i_radius,
3221 const unsigned j_radius,
3222 vil_image_view<float> const& intensity,
3223 const float range, const unsigned bins)
3224 {
3225 bsta_histogram<float> hi(range, bins);
3226 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3227 for (int dj = -jr; dj<=jr; ++dj)
3228 for (int di = -ir; di<=ir; ++di)
3229 {
3230 float inten = intensity(i+di, j+dj);
3231 hi.upcount(inten, 1.0f);
3232 }
3233 return hi.entropy();
3234 }
3235
3236 //: Compute the entropy of the gradient direction of a region
3237 // Note no bounds checking!
entropy_g(const unsigned i,const unsigned j,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & gradx,vil_image_view<float> const & grady,const float range,const unsigned bins)3238 float brip_vil_float_ops::entropy_g(const unsigned i, const unsigned j,
3239 const unsigned i_radius,
3240 const unsigned j_radius,
3241 vil_image_view<float> const& gradx,
3242 vil_image_view<float> const& grady,
3243 const float range, const unsigned bins)
3244 {
3245 bsta_histogram<float> hg(range, bins);
3246 static const auto deg_rad = (float)(vnl_math::deg_per_rad);
3247 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3248 for (int dj = -jr; dj<=jr; ++dj)
3249 for (int di = -ir; di<=ir; ++di)
3250 {
3251 float Ix = gradx(i+di, j+dj), Iy = grady(i+di, j+dj);
3252 float ang = deg_rad*std::atan2(Iy, Ix) + 180.0f;
3253 float mag = std::abs(Ix)+std::abs(Iy);
3254 hg.upcount(ang, mag);
3255 }
3256 return hg.entropy();
3257 }
3258
3259 //: Compute the hue and saturation entropy of a region about the specified pixel
3260 // Note no bounds checking!
entropy_hs(const unsigned i,const unsigned j,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & hue,vil_image_view<float> const & sat,const float range,const unsigned bins)3261 float brip_vil_float_ops::entropy_hs(const unsigned i, const unsigned j,
3262 const unsigned i_radius,
3263 const unsigned j_radius,
3264 vil_image_view<float> const& hue,
3265 vil_image_view<float> const& sat,
3266 const float range, const unsigned bins)
3267 {
3268 bsta_histogram<float> hg(range, bins);
3269 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3270 for (int dj = -jr; dj<=jr; ++dj)
3271 for (int di = -ir; di<=ir; ++di)
3272 {
3273 float h = hue(i+di, j+dj), s = sat(i+di, j+dj);
3274 hg.upcount(h, s);
3275 }
3276 return hg.entropy();
3277 }
3278
3279 vil_image_view<float>
entropy(const unsigned i_radius,const unsigned j_radius,const unsigned step,vil_image_resource_sptr const & img,const float sigma,const unsigned bins,const bool intensity,const bool gradient,const bool ihs)3280 brip_vil_float_ops::entropy(const unsigned i_radius,
3281 const unsigned j_radius,
3282 const unsigned step,
3283 vil_image_resource_sptr const& img,
3284 const float sigma,
3285 const unsigned bins,
3286 const bool intensity,
3287 const bool gradient,
3288 const bool ihs)
3289 {
3290 vil_image_view<float> ent;
3291 if (!intensity&&!gradient&&!ihs)
3292 {
3293 std::cout << "In brip_vil_float_ops::entropy(.) - No computation to do\n";
3294 return ent;
3295 }
3296
3297 vil_image_view<float> fimage = brip_vil_float_ops::convert_to_float(img);
3298 vil_image_view<float> gimage =
3299 brip_vil_float_ops::gaussian(fimage, sigma);
3300
3301 unsigned ni = img->ni(), nj = img->nj();
3302 ent.set_size(ni/step+1, nj/step+1);
3303 ent.fill(0.0f);
3304 if (intensity)
3305 for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
3306 for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
3307 ent(i/step,j/step) =
3308 brip_vil_float_ops::entropy_i(i, j, i_radius, j_radius, gimage, 255.0f, bins);
3309
3310 if (gradient)
3311 {
3312 vil_image_view<float> grad_x, grad_y;
3313 grad_x.set_size(ni, nj);
3314 grad_y.set_size(ni, nj);
3315 brip_vil_float_ops::gradient_3x3 (gimage , grad_x , grad_y);
3316 for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
3317 for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
3318 ent(i/step,j/step) +=
3319 brip_vil_float_ops::entropy_g(i, j, i_radius, j_radius,
3320 grad_x, grad_y);
3321 }
3322 if (ihs&&img->nplanes()==3)
3323 {
3324 vil_image_view<float> inten, hue, sat;
3325 vil_image_view<vil_rgb<vxl_byte> > cimage = img->get_view();
3326 brip_vil_float_ops::convert_to_IHS(cimage, inten, hue, sat);
3327 for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
3328 for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
3329 ent(i/step,j/step) +=
3330 brip_vil_float_ops::entropy_hs(i, j, i_radius, j_radius,
3331 hue, sat);
3332 }
3333 return ent;
3334 }
3335
minfo_i(const unsigned i0,const unsigned j0,const unsigned i1,const unsigned j1,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & int0,vil_image_view<float> const & int1,const float range,const unsigned bins)3336 float brip_vil_float_ops::minfo_i(const unsigned i0, const unsigned j0,
3337 const unsigned i1, const unsigned j1,
3338 const unsigned i_radius,
3339 const unsigned j_radius,
3340 vil_image_view<float> const& int0,
3341 vil_image_view<float> const& int1,
3342 const float range, const unsigned bins)
3343 {
3344 bsta_histogram<float> hi0(range, bins);
3345 bsta_histogram<float> hi1(range, bins);
3346 bsta_joint_histogram<float> hji(range, bins);
3347 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3348 for (int dj = -jr; dj<=jr; ++dj)
3349 for (int di = -ir; di<=ir; ++di)
3350 {
3351 float inten0 = int0(i0+di, j0+dj);
3352 float inten1 = int1(i1+di, j1+dj);
3353 hi0.upcount(inten0, 1.0f);
3354 hi1.upcount(inten1, 1.0f);
3355 hji.upcount(inten0, 1.0f, inten1, 1.0f);
3356 }
3357 float H0 = hi0.entropy();
3358 float H1 = hi1.entropy();
3359 float HJ = hji.entropy();
3360 float minfo_i = H0 + H1 - HJ;
3361 #ifdef DEBUG
3362 if (minfo<0)
3363 std::cout << "intensity MI LT 0 " << minfo << std::endl;
3364 #endif
3365 return minfo_i;
3366 }
3367
minfo_g(const unsigned i0,const unsigned j0,const unsigned i1,const unsigned j1,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & gradx0,vil_image_view<float> const & grady0,vil_image_view<float> const & gradx1,vil_image_view<float> const & grady1,const float range,const unsigned bins)3368 float brip_vil_float_ops::minfo_g(const unsigned i0, const unsigned j0,
3369 const unsigned i1, const unsigned j1,
3370 const unsigned i_radius,
3371 const unsigned j_radius,
3372 vil_image_view<float> const& gradx0,
3373 vil_image_view<float> const& grady0,
3374 vil_image_view<float> const& gradx1,
3375 vil_image_view<float> const& grady1,
3376 const float range, const unsigned bins)
3377 {
3378 bsta_histogram<float> hg0(range, bins);
3379 bsta_histogram<float> hg1(range, bins);
3380 bsta_joint_histogram<float> hjg(range, bins);
3381 static const auto deg_rad = (float)(vnl_math::deg_per_rad);
3382 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3383 for (int dj = -jr; dj<=jr; ++dj)
3384 for (int di = -ir; di<=ir; ++di)
3385 {
3386 float Ix0 = gradx0(i0+di, j0+dj), Iy0 = grady0(i0+di, j0+dj);
3387 float Ix1 = gradx1(i1+di, j1+dj), Iy1 = grady1(i1+di, j1+dj);
3388 float ang0 = deg_rad*std::atan2(Iy0, Ix0) + 180.0f;
3389 float ang1 = deg_rad*std::atan2(Iy1, Ix1) + 180.0f;
3390 float mag0 = std::abs(Ix0)+std::abs(Iy0);
3391 float mag1 = std::abs(Ix1)+std::abs(Iy1);
3392 hg0.upcount(ang0, mag0);
3393 hg1.upcount(ang1, mag1);
3394 hjg.upcount(ang0, mag0, ang1, mag1);
3395 }
3396 float H0 = hg0.entropy();
3397 float H1 = hg1.entropy();
3398 float HJ = hjg.entropy();
3399 float minfo_g = H0 + H1 - HJ;
3400 #ifdef DEBUG
3401 if (minfo<0)
3402 std::cout << "gradient MI LT 0 " << minfo << std::endl;
3403 #endif
3404 return minfo_g;
3405 }
3406
minfo_hs(const unsigned i0,const unsigned j0,const unsigned i1,const unsigned j1,const unsigned i_radius,const unsigned j_radius,vil_image_view<float> const & hue0,vil_image_view<float> const & sat0,vil_image_view<float> const & hue1,vil_image_view<float> const & sat1,const float range,const unsigned bins)3407 float brip_vil_float_ops::minfo_hs(const unsigned i0, const unsigned j0,
3408 const unsigned i1, const unsigned j1,
3409 const unsigned i_radius,
3410 const unsigned j_radius,
3411 vil_image_view<float> const& hue0,
3412 vil_image_view<float> const& sat0,
3413 vil_image_view<float> const& hue1,
3414 vil_image_view<float> const& sat1,
3415 const float range, const unsigned bins)
3416 {
3417 bsta_histogram<float> hh0(range, bins);
3418 bsta_histogram<float> hh1(range, bins);
3419 bsta_joint_histogram<float> hjh(range, bins);
3420 int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
3421 for (int dj = -jr; dj<=jr; ++dj)
3422 for (int di = -ir; di<=ir; ++di)
3423 {
3424 float h0 = hue0(i0+di, j0+dj), s0 = sat0(i0+di, j0+dj);
3425 float h1 = hue1(i1+di, j1+dj), s1 = sat1(i1+di, j1+dj);
3426 hh0.upcount(h0, s0);
3427 hh1.upcount(h1, s1);
3428 hjh.upcount(h0, s0, h1, s1);
3429 }
3430 float H0 = hh0.entropy();
3431 float H1 = hh1.entropy();
3432 float HJ = hjh.entropy();
3433 float minfo_h = H0 + H1 - HJ;
3434 #ifdef DEBUG
3435 if (minfo<0)
3436 std::cout << "color MI LT 0 " << minfo << std::endl;
3437 #endif
3438 return minfo_h;
3439 }
3440
minfo(const unsigned i_radius,const unsigned j_radius,const unsigned step,vil_image_resource_sptr const & img0,vil_image_resource_sptr const & img1,vil_image_view<float> & MI0,vil_image_view<float> & MI1,const float sigma,const bool intensity,const bool gradient,const bool ihs)3441 bool brip_vil_float_ops::minfo(const unsigned i_radius,
3442 const unsigned j_radius,
3443 const unsigned step,
3444 vil_image_resource_sptr const& img0,
3445 vil_image_resource_sptr const& img1,
3446 vil_image_view<float>& MI0,
3447 vil_image_view<float>& MI1,
3448 const float sigma,
3449 const bool intensity,
3450 const bool gradient,
3451 const bool ihs)
3452 {
3453 if (!intensity&&!gradient&&!ihs)
3454 {
3455 std::cout << "In brip_vil_float_ops::minforopy(.) - No computation to do\n";
3456 return false;
3457 }
3458
3459 vil_image_view<float> fimage0 = brip_vil_float_ops::convert_to_float(img0);
3460 vil_image_view<float> gimage0 =
3461 brip_vil_float_ops::gaussian(fimage0, sigma);
3462
3463 vil_image_view<float> fimage1 = brip_vil_float_ops::convert_to_float(img1);
3464 vil_image_view<float> gimage1 =
3465 brip_vil_float_ops::gaussian(fimage1, sigma);
3466
3467 unsigned ni0 = img0->ni(), nj0 = img0->nj();
3468 unsigned ni1 = img1->ni(), nj1 = img1->nj();
3469 unsigned ilimit = 2*i_radius +1, jlimit = 2*j_radius+1;
3470 if (ni0<ilimit||nj0<jlimit||ni1<ilimit||nj1<jlimit)
3471 {
3472 std::cout << "In brip_vil_float_ops::minfo(...) - image too small\n";
3473 return false;
3474 }
3475 MI0.set_size(ni0/step+1, nj0/step+1); MI0.fill(0.0f);
3476 MI1.set_size(ni1/step+1, nj1/step+1); MI1.fill(0.0f);
3477 if (intensity)
3478 for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
3479 for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
3480 for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
3481 for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
3482 {
3483 float minfo = brip_vil_float_ops::minfo_i(i0, j0,i1, j1,
3484 i_radius, j_radius,
3485 gimage0, gimage1);
3486 MI0(i0/step,j0/step) = minfo;
3487 MI1(i1/step,j1/step) = minfo;
3488 }
3489 if (gradient)
3490 {
3491 vil_image_view<float> grad_x0, grad_y0, grad_x1, grad_y1;
3492 grad_x0.set_size(ni0, nj0);
3493 grad_y0.set_size(ni0, nj0);
3494 grad_x1.set_size(ni1, nj1);
3495 grad_y1.set_size(ni1, nj1);
3496 brip_vil_float_ops::gradient_3x3 (gimage0 , grad_x0 , grad_y0);
3497 brip_vil_float_ops::gradient_3x3 (gimage1 , grad_x1 , grad_y1);
3498 for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
3499 for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
3500 for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
3501 for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
3502 {
3503 float minfo = brip_vil_float_ops::minfo_g(i0, j0,i1, j1,
3504 i_radius, j_radius,
3505 grad_x0, grad_y0,
3506 grad_x1, grad_y1);
3507 MI0(i0/step,j0/step) += minfo;
3508 MI1(i1/step,j1/step) += minfo;
3509 }
3510 }
3511 if (ihs&&img0->nplanes()==3&&img1->nplanes()==3)
3512 {
3513 vil_image_view<float> inten0, hue0, sat0;
3514 vil_image_view<float> inten1, hue1, sat1;
3515 vil_image_view<vil_rgb<vxl_byte> > cimage0 = img0->get_view();
3516 vil_image_view<vil_rgb<vxl_byte> > cimage1 = img1->get_view();
3517 brip_vil_float_ops::convert_to_IHS(cimage0, inten0, hue0, sat0);
3518 brip_vil_float_ops::convert_to_IHS(cimage1, inten1, hue1, sat1);
3519
3520 for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
3521 for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
3522 for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
3523 for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
3524 {
3525 float minfo = brip_vil_float_ops::minfo_hs(i0, j0,i1, j1,
3526 i_radius, j_radius,
3527 hue0, sat0,
3528 hue1, sat1);
3529 MI0(i0/step,j0/step) += minfo;
3530 MI1(i1/step,j1/step) += minfo;
3531 }
3532 }
3533 return true;
3534 }
3535
3536 // compute the average of the image intensity within the specified region
3537 float brip_vil_float_ops::
average_in_box(vil_image_view<float> const & v,vgl_box_2d<double> const & box)3538 average_in_box(vil_image_view<float> const& v, vgl_box_2d<double> const& box)
3539 {
3540 vgl_point_2d<double> p0 = box.min_point();
3541 auto i0 = static_cast<unsigned>(p0.x()), j0 = static_cast<unsigned>(p0.y());
3542 vgl_point_2d<double> p1 = box.max_point();
3543 auto i1 = static_cast<unsigned>(p1.x()), j1 = static_cast<unsigned>(p1.y());
3544 float n = 0.0f;
3545 float sum = 0.0f;
3546 for (unsigned i = i0; i<=i1; ++i)
3547 for (unsigned j = j0; j<=j1; ++j, ++n)
3548 sum += v(i,j);
3549 if (n>0)
3550 sum /= n;
3551 return sum;
3552 }
3553
3554 #if 0 // For now remove dependency on vimt. Save for illustration
3555 bool brip_vil_float_ops::vimt_homography(vil_image_view<float> const& curr_view,
3556 vgl_h_matrix_2d<double>const& H,
3557 vil_image_view<float>& output)
3558 {
3559 int bimg_ni;
3560 int bimg_nj;
3561
3562 int offset_i;
3563 int offset_j;
3564
3565 vbl_bounding_box<double,2> box;
3566
3567 unsigned ni = curr_view.ni(), nj = curr_view.nj();
3568 vimt_transform_2d p;
3569 vnl_double_3x3 Mh = H.get_matrix();
3570 vnl_double_2x3 M;
3571 for (unsigned c = 0; c<3; ++c)
3572 for (unsigned r = 0; r<2; ++r)
3573 M[r][c] = Mh[r][c]/Mh[2][2];
3574 p.set_affine(M);
3575 box.update(p(0,0).x(),p(0,0).y());
3576 box.update(p(0,nj).x(),p(0,nj).y());
3577 box.update(p(ni,0).x(),p(ni,0).y());
3578 box.update(p(ni,nj).x(),p(ni,nj).y());
3579
3580 bimg_ni=(int)std::ceil(box.max()[0]-box.min()[0]);
3581 bimg_nj=(int)std::ceil(box.max()[1]-box.min()[1]);
3582
3583 offset_i=(int)std::ceil(0-box.min()[0]);
3584 offset_j=(int)std::ceil(0-box.min()[1]);
3585
3586 vimt_image_2d_of<float> sample_im;
3587 vgl_point_2d<double> q(-offset_i,-offset_j);
3588 vgl_vector_2d<double> u(1,0);
3589 vgl_vector_2d<double> v(0,1);
3590
3591 vimt_image_2d_of<float> curr_img(curr_view,p);
3592 vimt_resample_bilin(curr_img,sample_im,q,u,v,bimg_ni,bimg_nj);
3593 output = sample_im.image();
3594 return true;
3595 }
3596
3597 #endif // 0
3598
scan_region(const vil_image_resource_sptr & img,const vgl_polygon<double> & poly,float & min,float & max)3599 std::vector<float> brip_vil_float_ops::scan_region(const vil_image_resource_sptr& img,
3600 const vgl_polygon<double>& poly,
3601 float& min,
3602 float& max)
3603 {
3604 std::vector<float> pixels;
3605 min = vnl_numeric_traits<float>::maxval;
3606 max = 0;
3607 unsigned np = img->nplanes();
3608 vgl_polygon_scan_iterator<double> si(poly, false);
3609 if (img->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
3610 {
3611 if (np==1) // single plane
3612 {
3613 for (si.reset(); si.next();)
3614 {
3615 auto j = static_cast<unsigned>(si.scany());
3616 for (int x = si.startx(); x<=si.endx(); ++x)
3617 {
3618 auto i = static_cast<unsigned>(x);
3619 vil_image_view<unsigned char> v = img->get_view(i, 1,j, 1);
3620 auto fv = static_cast<float>(v(0,0));
3621 if (fv<min) min = fv;
3622 if (fv>max) max = fv;
3623 pixels.push_back(fv);
3624 }
3625 }
3626 return pixels;
3627 }
3628 else
3629 {
3630 for (si.reset(); si.next();)
3631 {
3632 auto j = static_cast<unsigned>(si.scany());
3633 for (int x = si.startx(); x<=si.endx(); ++x)
3634 {
3635 auto i = static_cast<unsigned>(x);
3636 vil_image_view<unsigned char> v = img->get_view(i, 1,j, 1);
3637 float fv = 0;
3638 for (unsigned p = 0; p<np; ++p)
3639 fv += v(0,0,p);
3640 fv/=3;
3641 if (fv<min) min = fv;
3642 if (fv>max) max = fv;
3643 pixels.push_back(fv);
3644 }
3645 }
3646 }
3647 return pixels;
3648 }
3649 else if (img->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
3650 {
3651 if (np) // single plane
3652 {
3653 for (si.reset(); si.next();)
3654 {
3655 auto j = static_cast<unsigned>(si.scany());
3656 for (int x = si.startx(); x<=si.endx(); ++x)
3657 {
3658 auto i = static_cast<unsigned>(x);
3659 vil_image_view<unsigned short> v = img->get_view(i, 1,j, 1);
3660 auto fv = static_cast<float>(v(0,0));
3661 if (fv<min) min = fv;
3662 if (fv>max) max = fv;
3663 pixels.push_back(fv);
3664 }
3665 }
3666 return pixels;
3667 }
3668 else
3669 {
3670 for (si.reset(); si.next();)
3671 {
3672 auto j = static_cast<unsigned>(si.scany());
3673 for (int x = si.startx(); x<=si.endx(); ++x)
3674 {
3675 auto i = static_cast<unsigned>(x);
3676 vil_image_view<unsigned short> v = img->get_view(i, 1,j, 1);
3677 float fv = 0;
3678 for (unsigned p = 0; p<np; ++p)
3679 fv += v(0,0,p);
3680 fv/=3;
3681 if (fv<min) min = fv;
3682 if (fv>max) max = fv;
3683 pixels.push_back(fv);
3684 }
3685 }
3686 return pixels;
3687 }
3688 }
3689 std::cerr << "In brip_vil_float_ops::scan_region() - unknown format\n";
3690 return pixels;
3691 }
3692
3693 // It has been observed that color order is somewhat invariant to illumination
3694 // Compute an index based on color order. The tolerance determines if two
3695 // color bands are too close to determine order, i.e. they should be considered
3696 // equal instead
3697 // the two relations being considered are <, > and =, so the relationship
3698 // graph looks like:
3699 // \verbatim
3700 // G
3701 // / \.
3702 // > < = > < =
3703 // / \.
3704 // R - > < = - B
3705 // \endverbatim
3706 // Thus, there are three graph edges with each of three possible labels or
3707 // 9 possible order codes. An easy coding scheme is to use the top 6 bits of
3708 // the byte output pixel. The relationship is encoded as states of bit pairs
3709 // \verbatim
3710 // Color relations R*G R*B G*B * indicates > < = (1,2,3)
3711 // Bit indices 7,6 5,4 3,2
3712 // \endverbatim
3713 //
3714 vil_image_view<vxl_byte> brip_vil_float_ops::
color_order(vil_image_view<float> const & color_image,float eq_tol)3715 color_order(vil_image_view<float> const& color_image, float eq_tol)
3716 {
3717 unsigned ni = color_image.ni(), nj = color_image.nj();
3718 unsigned np = color_image.nplanes();
3719 vil_image_view<vxl_byte> temp;
3720 if (np!=3)
3721 return temp; // improper call
3722 temp.set_size(ni, nj);
3723 vxl_byte rg_eq=192, rg_lt=128, rg_gt=64;
3724 vxl_byte rb_eq=48, rb_lt=32, rb_gt=16;
3725 vxl_byte gb_eq=12, gb_lt=8, gb_gt=4;
3726 for (unsigned j = 0; j<nj; ++j)
3727 for (unsigned i = 0; i<ni; ++i) {
3728 float r=color_image(i,j,0), g=color_image(i,j,1), b=color_image(i,j,2);
3729 vxl_byte rg, rb, gb;
3730 // rg code
3731 if (std::fabs(r-g)<eq_tol)
3732 rg = rg_eq;
3733 else if (r<g)
3734 rg = rg_lt;
3735 else rg = rg_gt;
3736 // rb code
3737 if (std::fabs(r-b)<eq_tol)
3738 rb = rb_eq;
3739 else if (r<b)
3740 rb = rb_lt;
3741 else rb = rb_gt;
3742 // gb code
3743 if (std::fabs(g-b)<eq_tol)
3744 gb = gb_eq;
3745 else if (g<b)
3746 gb = gb_lt;
3747 else gb = gb_gt;
3748 auto v = static_cast<vxl_byte>(rg|rb|gb); // bitwise or
3749 temp(i,j) = v;
3750 }
3751 return temp;
3752 }
3753
3754 #if 0 // only used in commented-out part of brip_vil_float_ops::extrema()
3755 static double zs(double x)
3756 { if (x < 0.0001 && x > -0.0001) return 0; else return x; }
3757 #endif // 0
3758
brip_vil_rot_gauss(double x,double y,double sigma_x,double sigma_y,double theta)3759 static double brip_vil_rot_gauss(double x, double y,
3760 double sigma_x, double sigma_y, double theta)
3761 {
3762 double theta_rad = theta*vnl_math::pi_over_180;
3763 double s = std::sin(theta_rad), c = std::cos(theta_rad);
3764 double sxr = 1.0/sigma_x, syr = 1.0/sigma_y;
3765 double ax = (c*x + s*y)*sxr, ay = (-s*x + c*y)*syr;
3766 double ret = std::exp(-0.5*(ax*ax + ay*ay));
3767 return ret*sxr*syr/vnl_math::twopi;
3768 }
3769
3770 // revert angle to the range [-90, 90]
extrema_revert_angle(float angle)3771 float brip_vil_float_ops::extrema_revert_angle(float angle)
3772 {
3773 if (angle>=0.0f) {
3774 if (angle>90.0f&&angle<=270.0f)
3775 angle -= 180.0f;
3776 else if (angle>270.0f&&angle<=360.0f)
3777 angle -= 360.0f;
3778 }
3779 if (angle<0.0f) {
3780 if (angle<-90.0f&&angle>=-270.0f)
3781 angle += 180.0f;
3782 else if (angle<-270.0f&&angle>=-360.0f)
3783 angle += 360.0f;
3784 }
3785 return angle;
3786 }
3787
3788 // if scale_invariant = true then the response mag is independent of lambda1
extrema_kernel_mask(float lambda0,float lambda1,float theta,vbl_array_2d<float> & kernel,vbl_array_2d<bool> & mask,float cutoff_percentage,bool scale_invariant)3789 void brip_vil_float_ops::extrema_kernel_mask(float lambda0, float lambda1,
3790 float theta,
3791 vbl_array_2d<float>& kernel,
3792 vbl_array_2d<bool>& mask, float cutoff_percentage, bool scale_invariant)
3793 {
3794 theta = extrema_revert_angle(theta); // in range [-90, 90]
3795 // convert theta to radians
3796 double theta_rad = theta*vnl_math::pi_over_180;
3797 double s = std::sin(theta_rad), c = std::cos(theta_rad);
3798 double s0 = lambda0, s1 = lambda1;
3799 double s1sq = 1.0/(s1*s1);
3800 // determine the size of the array
3801 double max_v = brip_vil_rot_gauss(0, 0, s0, s1, 0);
3802 double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
3803 unsigned ru = 0;
3804 while (brip_vil_rot_gauss(ru, 0, s0, s1, 0) >= cutoff)
3805 ++ru;
3806 unsigned rv = 0;
3807 while (brip_vil_rot_gauss(0, rv, s0, s1, 0) >= cutoff)
3808 ++rv;
3809
3810 // rotate to get bounds
3811 int ri = static_cast<int>(std::fabs(ru*c+rv*s) +0.5);
3812 int rj = static_cast<int>(std::fabs(ru*s+rv*c) +0.5);
3813 if (s<0) {
3814 ri = static_cast<int>(std::fabs(ru*c-rv*s) +0.5);
3815 rj = static_cast<int>(std::fabs(ru*s-rv*c) +0.5);
3816 }
3817
3818 unsigned ncols = 2*ri +1, nrows = 2*rj+1;
3819 vbl_array_2d<double> coef(nrows, ncols);
3820 mask.resize(nrows,ncols);
3821 double residual = 0.0, total = 0.0;
3822 for (int ry = -rj; ry<=rj; ++ry)
3823 for (int rx = -ri; rx<=ri; ++rx)
3824 {
3825 double g = brip_vil_rot_gauss(rx, ry, s0, s1, theta);
3826 mask[ry+rj][rx+ri] = g>=cutoff;
3827 double temp = (-s*rx + c*ry);
3828 temp = temp*temp*s1sq;
3829 double v = (temp -1)*g*(scale_invariant?1.0:s1sq);
3830 if (g<cutoff) v = 0.0;
3831 coef[ry+rj][rx+ri] = v;
3832 residual+=v;
3833 total += std::fabs(v);
3834 }
3835 double cor = 0.0;
3836 if (total != 0)
3837 cor = -residual/total;
3838 kernel.resize(nrows, ncols);
3839 // correct any residual offset in coefficients
3840 // distribute proportionally to coefficient magnitude
3841 for (unsigned j = 0; j<nrows; ++j)
3842 for (unsigned i = 0; i<ncols; ++i)
3843 {
3844 double v = std::fabs(coef[j][i]);
3845 coef[j][i]+=v*cor;
3846 kernel[j][i]=static_cast<float>(coef[j][i]);
3847 }
3848 }
3849
3850 //:
3851 // \p theta must be given in degrees.
3852 // Scale invariant means that the response is independent of the
3853 // \a sigma_y of the unrotated operator, i.e. the direction of the derivative
gaussian_kernel_mask(float lambda,vbl_array_2d<float> & kernel,vbl_array_2d<bool> & mask,float cutoff_percentage)3854 void brip_vil_float_ops::gaussian_kernel_mask(float lambda, vbl_array_2d<float>& kernel, vbl_array_2d<bool>& mask, float cutoff_percentage)
3855 {
3856 double s1 = lambda;
3857 double s1sq = 1.0/(s1*s1);
3858 // determine the size of the array
3859 double max_v = brip_vil_rot_gauss(0, 0, s1, s1, 0);
3860 double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
3861 unsigned r = 0;
3862 while (brip_vil_rot_gauss(r, 0, s1, s1, 0) >= cutoff)
3863 ++r;
3864
3865 int ri = static_cast<int>(std::fabs((float)r) +0.5);
3866
3867 unsigned n = 2*ri +1;
3868 mask.resize(n,n);
3869 kernel.resize(n, n);
3870 for (int ry = -ri; ry<=ri; ++ry)
3871 for (int rx = -ri; rx<=ri; ++rx)
3872 {
3873 double g = brip_vil_rot_gauss(rx, ry, s1, s1, 0);
3874 mask[ry+ri][rx+ri] = g>=cutoff;
3875 double temp = ry*s1sq;
3876 double v = (temp -1)*g;
3877 if (g<cutoff) v = 0.0;
3878 kernel[ry+ri][rx+ri] = static_cast<float>(std::fabs(v));
3879 }
3880 }
3881
3882 //: return a square mask and the coefficient matrix [kernel] for a symmetric gaussian distribution
gaussian_kernel_square_mask(float lambda,vbl_array_2d<float> & kernel,vbl_array_2d<bool> & mask,float cutoff_percentage)3883 void brip_vil_float_ops::gaussian_kernel_square_mask(float lambda, vbl_array_2d<float>& kernel, vbl_array_2d<bool>& mask, float cutoff_percentage)
3884 {
3885 double s1 = lambda;
3886 double s1sq = 1.0/(s1*s1);
3887 // determine the size of the array
3888 double max_v = brip_vil_rot_gauss(0, 0, s1, s1, 0);
3889 double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
3890 unsigned r = 0;
3891 while (brip_vil_rot_gauss(r, 0, s1, s1, 0) >= cutoff)
3892 ++r;
3893
3894 int ri = static_cast<int>(std::fabs((float)r) +0.5);
3895
3896 unsigned n = 2*ri +1;
3897 vbl_array_2d<double> coef(n, n);
3898 mask.resize(n,n);
3899 kernel.resize(n, n);
3900 for (int ry = -ri; ry<=ri; ++ry)
3901 for (int rx = -ri; rx<=ri; ++rx)
3902 {
3903 double g = brip_vil_rot_gauss(rx, ry, s1, s1, 0);
3904 mask[ry+ri][rx+ri] = true;
3905 double temp = ry*s1sq;
3906 kernel[ry+ri][rx+ri] = static_cast<float>(std::fabs((temp -1)*g));
3907 }
3908 }
3909
3910
3911 // Compute the standard deviation of an operator response
3912 // given the image intensity standard deviation at each pixel
3913 vil_image_view<float> brip_vil_float_ops::
std_dev_operator(vil_image_view<float> const & sd_image,vbl_array_2d<float> const & kernel)3914 std_dev_operator(vil_image_view<float> const& sd_image,
3915 vbl_array_2d<float> const& kernel)
3916 {
3917 unsigned ni = sd_image.ni(), nj = sd_image.nj();
3918 unsigned nrows = kernel.rows(), ncols = kernel.cols();
3919 int rrad = (nrows-1)/2, crad = (ncols-1)/2;
3920 float sd_min, sd_max;
3921 vil_math_value_range(sd_image, sd_min, sd_max);
3922
3923 vil_image_view<float> res(ni, nj);
3924 res.fill(sd_max);
3925
3926 vbl_array_2d<float>ksq(nrows, ncols);
3927 for (unsigned r = 0; r<nrows; ++r)
3928 for (unsigned c = 0; c<ncols; ++c)
3929 ksq[r][c] = kernel[r][c]*kernel[r][c];
3930
3931 for (int j = rrad; j<static_cast<int>(nj-rrad); ++j)
3932 for (int i = crad; i<static_cast<int>(ni-crad); ++i)
3933 {
3934 float sum = 0;
3935 for (int jj = -rrad; jj<=rrad; ++jj)
3936 for (int ii = -crad; ii<=crad; ++ii) {
3937 float sd_sq = sd_image(i+ii, j+jj);
3938 sd_sq *= sd_sq;
3939 sum += sd_sq*ksq[jj+rrad][ii+crad];
3940 }
3941 res(i,j) = std::sqrt(sum);
3942 }
3943 return res;
3944 }
3945
3946 // Compute the standard deviation of an operator response
3947 // given the image intensity standard deviation at each pixel
3948 // uses a modified formula to compute std_dev
3949 vil_image_view<float> brip_vil_float_ops::
std_dev_operator_method2(vil_image_view<float> const & sd_image,vbl_array_2d<float> const & kernel)3950 std_dev_operator_method2(vil_image_view<float> const& sd_image,
3951 vbl_array_2d<float> const& kernel)
3952 {
3953 unsigned ni = sd_image.ni(), nj = sd_image.nj();
3954 unsigned nrows = kernel.rows(), ncols = kernel.cols();
3955 int rrad = (nrows-1)/2, crad = (ncols-1)/2;
3956 float sd_min, sd_max;
3957 vil_math_value_range(sd_image, sd_min, sd_max);
3958
3959 vil_image_view<float> res(ni, nj);
3960 res.fill(sd_max);
3961
3962 for (int j = rrad; j<static_cast<int>(nj-rrad); ++j)
3963 for (int i = crad; i<static_cast<int>(ni-crad); ++i)
3964 {
3965 float sum = 0;
3966 for (int jj = -rrad; jj<=rrad; ++jj)
3967 for (int ii = -crad; ii<=crad; ++ii) {
3968 float sd = sd_image(i+ii, j+jj);
3969 sum += (float)(sd*std::abs(kernel[jj+rrad][ii+crad]));
3970 }
3971 res(i,j) = sum;
3972 }
3973 return res;
3974 }
3975
3976 vil_image_view<float>
extrema(vil_image_view<float> const & input,float lambda0,float lambda1,float theta,bool bright,bool mag_only,bool output_response_mask,bool signed_response,bool scale_invariant,bool non_max_suppress,float cutoff_per)3977 brip_vil_float_ops::extrema(vil_image_view<float> const& input,
3978 float lambda0, float lambda1, float theta,
3979 bool bright, bool mag_only,
3980 bool output_response_mask,
3981 bool signed_response, bool scale_invariant,
3982 bool non_max_suppress, float cutoff_per)
3983 {
3984 assert(!(mag_only&&signed_response));
3985 vbl_array_2d<float> fa;
3986 vbl_array_2d<bool> mask;
3987 brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta,
3988 fa, mask, cutoff_per,
3989 scale_invariant);
3990 unsigned nrows = fa.rows(), ncols = fa.cols();
3991 int rj = (nrows-1)/2, ri = (ncols-1)/2;
3992 vbl_array_2d<double> coef(nrows,ncols);
3993 for (unsigned r = 0; r<nrows; ++r)
3994 for (unsigned c = 0; c<ncols; ++c)
3995 coef[r][c]=fa[r][c];
3996
3997 #if 0
3998 // (unused) double max_v = brip_vil_rot_gauss(0, 0, lambda0, lambda1, 0);
3999 // (unused) double cutoff=static_cast<float>(max_v*0.01); // 1% tails removed
4000 std::cout << "\ngauss ={";
4001 for (unsigned j = 0; j<nrows; ++j) {
4002 std::cout << '{';
4003 for (unsigned i = 0; i<ncols-1; ++i)
4004 std::cout << zs(coef[j][i]) << ',';
4005 if (j != nrows-1)
4006 std::cout << zs(coef[j][ncols-1]) << "},";
4007 else
4008 std::cout << zs(coef[j][ncols-1]) << '}';
4009 }
4010
4011 std::cout << "};\n\nmask ={";
4012 for (unsigned j = 0; j<nrows; ++j) {
4013 std::cout << '{';
4014 for (unsigned i = 0; i<ncols-1; ++i)
4015 std::cout << mask[j][i] << ',';
4016 if (j != nrows-1)
4017 std::cout << mask[j][ncols-1] << "},";
4018 else
4019 std::cout << mask[j][ncols-1] << '}';
4020 }
4021 std::cout << "};" << std::flush;
4022 #endif
4023 // compute response image
4024 unsigned ni = input.ni(), nj = input.nj();
4025 vil_image_view<float> temp(ni, nj);
4026 vil_image_view<float> temp2(ni, nj);
4027 temp.fill(0.0f); temp2.fill(0.0f);
4028 for (unsigned j = rj; j<(nj-rj); j++)
4029 for (unsigned i = ri; i<(ni-ri); i++) {
4030 double sum = 0;
4031 for (int jj=-rj; jj<=rj; ++jj)
4032 for (int ii=-ri; ii<=ri; ++ii)
4033 if (mask[jj+rj][ii+ri])
4034 sum += coef[jj+rj][ii+ri]*input(i+ii, j+jj);
4035 temp2(i,j) = static_cast<float>(sum);
4036 if (mag_only) {
4037 temp(i,j) = static_cast<float>(std::fabs(sum));
4038 }
4039 else if (bright) { // coefficients are negative at center
4040 if (sum<0) temp(i,j) = static_cast<float>(-sum);
4041 }
4042 else {
4043 if (sum>0) temp(i,j) = static_cast<float>(sum);
4044 }
4045 }
4046 // non-max suppression
4047 vil_image_view<float> res(temp);
4048 if (non_max_suppress)
4049 for (unsigned j = rj; j<(nj-rj); j++)
4050 for (unsigned i = ri; i<(ni-ri); i++)
4051 {
4052 float cv = temp(i,j);
4053 for (int jj=-rj; jj<=rj; ++jj)
4054 for (int ii=-ri; ii<=ri; ++ii)
4055 if ((ii!=0||jj!=0) && mask[jj+rj][ii+ri] && temp(i+ii, j+jj)>cv)
4056 res(i,j)=0.0f;
4057 }
4058 if (!output_response_mask&&!signed_response)
4059 return res;
4060 unsigned np = 2;
4061 if (output_response_mask&&signed_response)
4062 np = 3;
4063 //response plane and mask plane and/or signed response
4064 vil_image_view<float> res_mask(ni, nj, np);
4065 res_mask.fill(0.0f);
4066 // always copy response to plane 0
4067 for (unsigned j = rj; j<(nj-rj); j++)
4068 for (unsigned i = ri; i<(ni-ri); i++)
4069 {
4070 float rv = res(i,j);
4071 res_mask(i,j,0)=rv;
4072 }
4073 // copy mask plane to plane 1 (or not)
4074 if (output_response_mask) {
4075 for (unsigned j = rj; j<(nj-rj); j++)
4076 for (unsigned i = ri; i<(ni-ri); i++)
4077 {
4078 float rv = res(i,j);
4079 if (!rv)
4080 continue;
4081 for (int jj=-rj; jj<=rj; ++jj)
4082 for (int ii=-ri; ii<=ri; ++ii)
4083 if (mask[jj+rj][ii+ri])
4084 if (rv>res_mask(i+ii,j+jj,1))
4085 res_mask(i+ii,j+jj,1) = rv;
4086 }
4087 }
4088 // copy signed response to plane 1 if no mask plane, otherwise plane 2
4089 if (signed_response) {
4090 unsigned p = np-1;
4091 for (unsigned j = rj; j<(nj-rj); j++)
4092 for (unsigned i = ri; i<(ni-ri); i++)
4093 res_mask(i,j,p) = temp2(i,j);
4094 }
4095 return res_mask;
4096 }
4097
4098
4099 //: Find anisotropic intensity extrema at a range of orientations and return the maximal response at the best orientation. Theta interval is in degrees
4100 // if lambda0 == lambda1 then reduces to the normal extrema operator
4101 vil_image_view<float> brip_vil_float_ops::
extrema_rotational(vil_image_view<float> const & input,float lambda0,float lambda1,float theta_interval,bool bright,bool mag_only,bool signed_response,bool scale_invariant,bool non_max_suppress,float cutoff_per)4102 extrema_rotational(vil_image_view<float> const& input, float lambda0,
4103 float lambda1, float theta_interval, bool bright,
4104 bool mag_only, bool signed_response,
4105 bool scale_invariant, bool non_max_suppress,
4106 float cutoff_per)
4107 {
4108 //invalid response
4109 assert(!(mag_only&&signed_response));
4110
4111 unsigned ni = input.ni();
4112 unsigned nj = input.nj();
4113
4114 vil_image_view<float> output(ni, nj, 3); // the first plane is the response strength, the second is the best angle and the third is the mask for that angle
4115 output.fill(0.0f);
4116
4117 if (lambda0 == lambda1) { // if isotropic reduces to normal extrema calculation (perfect orientation symmetry)
4118 vil_image_view<float> res_mask_current =
4119 extrema(input, lambda0, lambda1, 0.0f, bright, mag_only,
4120 true, signed_response, scale_invariant,
4121 non_max_suppress, cutoff_per);
4122 unsigned np = res_mask_current.nplanes();
4123 for (unsigned j = 0; j<nj; j++)
4124 for (unsigned i = 0; i<ni; i++) {
4125 if (!signed_response)
4126 output(i,j,0) = res_mask_current(i,j,0); // return the non-max suppressed for angle 0
4127 else
4128 output(i,j,0) = res_mask_current(i,j,np-1);
4129 output(i,j,1) = 0.0f;
4130 output(i,j,2) = res_mask_current(i,j,1);
4131 }
4132 return output;
4133 }
4134
4135 // The kernel generator does not treat the x and y axis symmetrically, the method works correctly only when lambda0 > lambda1
4136 // Theoretically one can always call this method by switching the lambdas but the caller of the method should make this switch if needed hence the assertion
4137 if (lambda0 < lambda1) {
4138 std::cout << "In brip_vil_float_ops::extrema_rotational() - ERROR! rotational extrema operator requires lambda0 to be larger than lambda1! switch the lambdas and use the output angle map accordingly!\n";
4139 throw 0;
4140 }
4141 if (theta_interval < std::numeric_limits<float>::epsilon()) {
4142 std::cout << "In brip_vil_float_ops::extrema_rotational() - ERROR! theta_interval needs to be larger than 0!\n";
4143 throw 0;
4144 }
4145
4146 vil_image_view<float> res_img(ni, nj);
4147 vil_image_view<int> res_angle(ni, nj);
4148 res_img.fill(0.0f); res_angle.fill(0);
4149
4150 std::vector<float> angles;
4151 angles.push_back(-1.0f);
4152 for (float theta = 0.0f; theta < 180.0f; theta += theta_interval) { angles.push_back(theta); }
4153
4154 std::vector<vbl_array_2d<bool> > mask_vect(angles.size(), vbl_array_2d<bool>());
4155 int max_rji = 0;
4156 // elliptical operator has 180 degree rotational symmetry, so only the angles in the range [0,180] matter
4157 float theta = 0.0f; unsigned theta_i = 1;
4158 for ( ; theta < 180; theta += theta_interval, theta_i++)
4159 {
4160 // compute the response
4161 vbl_array_2d<float> fa; vbl_array_2d<bool> mask;
4162 brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta,
4163 fa, mask, cutoff_per, scale_invariant);
4164 mask_vect[theta_i] = mask;
4165 unsigned nrows = fa.rows(), ncols = fa.cols();
4166 int rj = (nrows-1)/2, ri = (ncols-1)/2;
4167 if (rj > max_rji) max_rji = rj;
4168 if (ri > max_rji) max_rji = ri;
4169 for (unsigned j = rj; j<(nj-rj); j++) {
4170 for (unsigned i = ri; i<(ni-ri); i++) {
4171 float res = 0.0f;
4172 double sum = 0;
4173 for (int jj=-rj; jj<=rj; ++jj)
4174 for (int ii=-ri; ii<=ri; ++ii)
4175 if (mask[jj+rj][ii+ri]) {
4176 sum += double(fa[jj+rj][ii+ri])*input(i+ii, j+jj);
4177 }
4178 if (mag_only) {
4179 res = static_cast<float>(std::fabs(sum));
4180 }
4181 else if (signed_response) {
4182 res = static_cast<float>(sum);
4183 }
4184 else if (bright) { // coefficients are negative at center
4185 if (sum<0) res = static_cast<float>(-sum);
4186 }
4187 else {
4188 if (sum>0) res = static_cast<float>(sum);
4189 }
4190 // pick largest angle response magnitude
4191 if (std::fabs(res_img(i,j)) < std::fabs(res)) {
4192 res_img(i,j) = res;
4193 res_angle(i,j) = theta_i;
4194 output(i,j,0) = res;
4195 output(i,j,1) = theta_i;
4196 }
4197 }
4198 }
4199 std::cout << '.';
4200 }
4201 std::cout << '\n';
4202 //if (!non_max_suppress) return res_img;
4203 if (!non_max_suppress) return output;
4204 // now we have pixel-wise best angle, run the non-max suppression around each non-zero pixel using the angles mask
4205 vil_image_view<float> res(res_img);
4206
4207 for (unsigned j = max_rji; j<(nj-max_rji); j++)
4208 for (unsigned i = max_rji; i<(ni-max_rji); i++)
4209 {
4210 float cv = res_img(i,j);
4211 if (!cv)
4212 continue;
4213 int theta_i = res_angle(i,j);
4214 vbl_array_2d<bool> mask = mask_vect[theta_i];
4215 unsigned nrows = mask.rows(), ncols = mask.cols();
4216 int rj = (nrows-1)/2, ri = (ncols-1)/2;
4217
4218 bool max = true;
4219 for (int jj=-rj; jj<=rj; ++jj) {
4220 for (int ii=-ri; ii<=ri; ++ii) {
4221 if ((ii!=0||jj!=0) && mask[jj+rj][ii+ri] && res_img(i+ii, j+jj)>cv) {
4222 max = false;
4223 break;
4224 }
4225 }
4226 if (!max)
4227 break;
4228 }
4229 if (!max) {
4230 res(i,j) = 0.0f;
4231 res_angle(i, j) = 0; // the zeroth angle is -1.0f, so invalid
4232 }
4233 }
4234
4235 vil_image_view<float> res_mask(ni, nj);
4236 res_mask.fill(0.0f);
4237 // now all the non-zero elements in res are our responses, create the mask image using the angle info
4238 for (unsigned j = max_rji; j<(nj-max_rji); j++)
4239 for (unsigned i = max_rji; i<(ni-max_rji); i++)
4240 {
4241 float rv = res(i,j);
4242 if (!rv)
4243 continue;
4244 int theta_i = res_angle(i,j);
4245 // get the mask for this angle
4246 vbl_array_2d<bool> mask = mask_vect[theta_i];
4247 unsigned nrows = mask.rows(), ncols = mask.cols();
4248 int rj = (nrows-1)/2, ri = (ncols-1)/2;
4249
4250 for (int jj=-rj; jj<=rj; ++jj)
4251 for (int ii=-ri; ii<=ri; ++ii)
4252 if (mask[jj+rj][ii+ri])
4253 if (rv>res_mask(i+ii,j+jj))
4254 res_mask(i+ii,j+jj) = rv;
4255 }
4256
4257 // now prepare the output accordingly
4258 for (unsigned j = max_rji; j<(nj-max_rji); j++)
4259 for (unsigned i = max_rji; i<(ni-max_rji); i++)
4260 {
4261 output(i,j,0) = res(i,j);
4262 output(i,j,1) = angles[res_angle(i,j)];
4263 output(i,j,2) = res_mask(i,j);
4264 }
4265 return output;
4266 }
4267
4268
4269 //: theta and phi are in radians
elu(float phi,float lamda0,float lambda1,float theta)4270 float brip_vil_float_ops::elu(float phi, float lamda0,
4271 float lambda1, float theta)
4272 {
4273 double sp = std::sin(phi), cp = std::cos(phi);
4274 double st = std::sin(theta), ct = std::cos(theta);
4275 double temp = 3.0*lamda0*cp*ct-3.0*lambda1*sp*st;
4276 return static_cast<float>(temp);
4277 }
4278
4279 //: theta and phi are in radians
elv(float phi,float lamda0,float lambda1,float theta)4280 float brip_vil_float_ops::elv(float phi, float lamda0,
4281 float lambda1, float theta)
4282 {
4283 double sp = std::sin(phi), cp = std::cos(phi);
4284 double st = std::sin(theta), ct = std::cos(theta);
4285 double temp = 3.0*lamda0*st*cp + 3.0*lambda1*sp*ct;
4286 return static_cast<float>(temp);
4287 }
4288
4289 //: theta and phi are in radians
4290 void brip_vil_float_ops::
max_inscribed_rect(float lambda0,float lambda1,float theta,float & u_rect,float & v_rect)4291 max_inscribed_rect(float lambda0, float lambda1, float theta,
4292 float& u_rect, float& v_rect)
4293 {
4294 float sign = -1.0f;
4295 if (theta<0) sign = 1.0f;
4296 auto th_rad = static_cast<float>(theta*vnl_math::pi_over_180);
4297 float maxa = 0.0f, max_phi = 0.0f;
4298 for (float phi = -float(vnl_math::pi); phi<=float(vnl_math::pi); phi+=0.01f)
4299 {
4300 float a = (1.0f+brip_vil_float_ops::elu(phi, lambda0, lambda1, th_rad));
4301 a *= (1.0f+ sign*brip_vil_float_ops::elv(phi, lambda0, lambda1, th_rad));
4302 if (a>maxa)
4303 {
4304 maxa = a;
4305 max_phi = phi;
4306 }
4307 }
4308 u_rect = brip_vil_float_ops::elu(max_phi, lambda0, lambda1, th_rad);
4309 v_rect = brip_vil_float_ops::elv(max_phi, lambda0, lambda1, th_rad);
4310 v_rect = std::fabs(v_rect);
4311 }
4312
eu(float lambda0,float cutoff_per,std::vector<float> & kernel)4313 static void eu(float lambda0, float cutoff_per, std::vector<float>& kernel)
4314 {
4315 kernel.clear();
4316 double l_sqi = 1.0/(lambda0*lambda0);
4317 double norm = vnl_math::one_over_sqrt2pi;
4318 norm /= lambda0;
4319 int r = static_cast<int>(std::sqrt((-2.0*std::log(cutoff_per))/l_sqi)+0.5);
4320 for (int i = -r; i<=r; ++i)
4321 {
4322 double k = std::exp(-0.5*i*i*l_sqi);
4323 k*=norm;
4324 kernel.push_back(static_cast<float>(k));
4325 }
4326 }
4327
ev(float lambda1,float cutoff_per,bool scale_invariant,std::vector<float> & kernel)4328 static void ev(float lambda1, float cutoff_per, bool scale_invariant,
4329 std::vector<float>& kernel)
4330 {
4331 kernel.clear();
4332 double l1_sqi = 1/(lambda1*lambda1);
4333 int r = static_cast<int>(std::sqrt((-2.0*std::log(cutoff_per))/l1_sqi)+0.5);
4334 double norm = scale_invariant ? 1/lambda1 : l1_sqi/lambda1;
4335 norm *= vnl_math::one_over_sqrt2pi;
4336 for (int i = -r; i<=r; ++i)
4337 {
4338 double s = i*i*l1_sqi;
4339 double k = std::exp(-0.5*s);
4340 k *= norm*(s-1.0);
4341 kernel.push_back(static_cast<float>(k));
4342 }
4343 }
4344
convolve_u(vil_image_view<float> const & image,std::vector<float> const & kernel,bool initialize=true)4345 static vil_image_view<float> convolve_u(vil_image_view<float> const& image,
4346 std::vector<float> const& kernel,
4347 bool initialize = true)
4348 {
4349 int nk = kernel.size();
4350 if (nk<3)
4351 return image;
4352 int ni = image.ni(), nj = image.nj();
4353 int ru = (nk-1)/2;
4354 vil_image_view<float> Iu(ni, nj);
4355 if (initialize) Iu.fill(0.0f);
4356 for (int j = 0; j<nj; ++j)
4357 for (int i = ru; i<ni-ru; ++i) {
4358 float sum = 0.0f;
4359 for (int k = -ru; k<=ru; ++k)
4360 {
4361 float v = image(i+k, j);
4362 sum += kernel[k+ru]*v;
4363 }
4364 Iu(i,j) = sum;
4365 }
4366 return Iu;
4367 }
4368
convolve_v(vil_image_view<float> const & image,std::vector<float> const & kernel,bool initialize=true)4369 static vil_image_view<float> convolve_v(vil_image_view<float> const& image,
4370 std::vector<float> const& kernel,
4371 bool initialize = true)
4372 {
4373 int nk = kernel.size();
4374 if (nk<3)
4375 return image;
4376 int ni = image.ni(), nj = image.nj();
4377 int rv = (nk-1)/2;
4378 vil_image_view<float> Iv(ni, nj);
4379 if (initialize) Iv.fill(0.0f);
4380 for (int j = rv; j<nj-rv; ++j)
4381 for (int i = 0; i<ni; ++i) {
4382 float sum = 0.0f;
4383 for (int k = -rv; k<=rv; ++k)
4384 {
4385 float v = image(i, j+k);
4386 sum += kernel[k+rv]*v;
4387 }
4388 Iv(i,j) = sum;
4389 }
4390 return Iv;
4391 }
4392
4393 // this function tracks the image origin during rotation
4394 // the calculations are the same as used in the homography method.
rotation_offset(int ni,int nj,float theta_deg,int i,int j,int & ti,int & tj)4395 static void rotation_offset(int ni, int nj, float theta_deg,
4396 int i, int j, int& ti, int& tj)
4397 {
4398 ti= 0; tj = 0;
4399 double deg_to_rad = vnl_math::pi_over_180;
4400 double rang = deg_to_rad*theta_deg;
4401 vgl_h_matrix_2d<double> H;
4402 H.set_identity().set_rotation(rang);
4403 vsol_box_2d_sptr input_roi;
4404 vsol_polygon_2d_sptr input_poly, output_poly;
4405 input_roi = new vsol_box_2d();
4406 input_roi->add_point(0, 0);
4407 input_roi->add_point(ni, nj);
4408 input_poly = bsol_algs::poly_from_box(input_roi);
4409 if (!bsol_algs::homography(input_poly, H, output_poly))
4410 return;
4411 vsol_box_2d_sptr temp = output_poly->get_bounding_box();
4412 vgl_h_matrix_2d<double> t; t.set_identity().set_translation(-temp->get_min_x(),-temp->get_min_y());
4413 vgl_homg_point_2d<double> torg = (t*H)*vgl_homg_point_2d<double>(double(i),double(j));
4414 ti = static_cast<int>(torg.x());
4415 tj = static_cast<int>(torg.y());
4416 }
4417
4418 vil_image_view<float> brip_vil_float_ops::
fast_extrema(vil_image_view<float> const & input,float lambda0,float lambda1,float theta,bool bright,bool mag_only,bool output_response_mask,bool signed_response,bool scale_invariant,bool non_max_suppress,float cutoff)4419 fast_extrema(vil_image_view<float> const& input, float lambda0, float lambda1,
4420 float theta, bool bright, bool mag_only,bool output_response_mask,
4421 bool signed_response, bool scale_invariant, bool non_max_suppress,
4422 float cutoff)
4423 {
4424 // invalid choice
4425 assert(!(mag_only&&signed_response));
4426 // the kernels for u and v
4427 std::vector<float> euk, evk;
4428 eu(lambda0, cutoff, euk);
4429 ev(lambda1, cutoff, scale_invariant, evk);
4430 vil_image_view<float> pre_res;
4431 int tuf=0, tvf=0, tui=0, tvi=0;
4432 vil_image_view<float> rot = input;
4433 bool ang90 = theta==90.0f||theta == 270.0f||
4434 theta ==-90.0f||theta == -270.0f;
4435 // rotate the image to the specified angle (minus due to flipped v axis)
4436 if (theta!=0&&!ang90)
4437 rot = brip_vil_float_ops::rotate(input, -theta);
4438 if (!ang90) {
4439 // convolve kernel across rows
4440 vil_image_view<float> rotu = convolve_u(rot, euk);
4441 // then convolve along columns
4442 vil_image_view<float> rotuv = convolve_v(rotu, evk);
4443 // rotate back
4444 if (theta!=0)
4445 {
4446 pre_res = brip_vil_float_ops::rotate(rotuv, +theta);
4447 // there is a border around the image to contain the rotated version
4448 // which should be removed.
4449 rotation_offset(input.ni(), input.nj(), -theta, 0, 0, tuf, tvf);
4450 rotation_offset(rotuv.ni(), rotuv.nj(), +theta, tuf, tvf, tui, tvi);
4451 }
4452 else pre_res = rotuv;
4453 }
4454 else {
4455 // convolve along columns
4456 vil_image_view<float> rotv = convolve_v(rot, euk);
4457 // convolve across rows
4458 vil_image_view<float> rotvu = convolve_u(rotv, evk);
4459 pre_res = rotvu;
4460 }
4461 int ni = input.ni(), nj = input.nj();
4462 vil_image_view<float> res(ni, nj);
4463 res.fill(0.0f);
4464 vil_image_view<float> temp2(ni, nj);
4465 temp2.fill(0.0f);
4466 for (int j = 0; j<nj; ++j)
4467 for (int i = 0; i<ni; ++i)
4468 {
4469 float v = pre_res(i + tui, j+ tvi);
4470 if (mag_only) {
4471 res(i,j) = std::fabs(v);
4472 }
4473 else if (bright) { // coefficients are negative at center
4474 if (v<0) res(i,j) = -v;
4475 }
4476 else {
4477 if (v>0) res(i,j) = v;
4478 }
4479 temp2(i,j) = v;// signed output
4480 }
4481 if (!non_max_suppress && !output_response_mask) {
4482 return signed_response ? temp2 : res;
4483 }
4484 // TO DO! handle non_max_suppress==false but response_mask==true
4485 // December 12, 2011 --- JLM
4486 //
4487 // the second derivative convolution is complete
4488 // non-maximal suppression
4489 // determine the inscribed rect in the response mask with max area
4490 float u_rect, v_rect;
4491 brip_vil_float_ops::max_inscribed_rect(lambda0,lambda1,theta,u_rect,v_rect);
4492 int ni_rect = static_cast<int>(u_rect);
4493 int nj_rect = static_cast<int>(v_rect);
4494 vil_image_view<float> res_sup(ni, nj);
4495 res_sup.fill(0.0f);
4496 for ( int i= ni_rect; i< ni-ni_rect; i+= ni_rect+1 )
4497 for ( int j = nj_rect; j < nj-nj_rect; j += nj_rect+1 ) {
4498 int ci= i, cj = j;
4499 bool find_response = false;
4500 for ( int di= 0; di<= ni_rect; di++ )
4501 for ( int dj = 0; dj <= nj_rect; dj++ )
4502 {
4503 float v = res(ci,cj);
4504 float dv = res(i+di,j+dj);
4505 if (v==0.0f&&dv==0.0f) continue;
4506 find_response = true;
4507 if ( v < dv )
4508 {
4509 ci= i+di;
4510 cj = j+dj;
4511 }
4512 }
4513 if (!find_response)
4514 continue;
4515 res_sup(ci,cj) = res(ci,cj);
4516 }
4517
4518 // determine the outside bounding box
4519 vbl_array_2d<float> kern;
4520 vbl_array_2d<bool> mask;
4521 brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta, kern, mask);
4522 int ri = (kern.cols()-1)/2, rj = (kern.rows()-1)/2;
4523 // go through the array again and search about each retained
4524 // local maximum using the outside bounding box
4525 // there will be at most one non-zero response to test in the inscribed rect
4526 for (int j = rj; j<(nj-rj); j++)
4527 for (int i = ri; i<(ni-ri); i++)
4528 {
4529 float v = res_sup(i,j);
4530 if (v==0.0f) continue;
4531 for (int jj=-rj; jj<=rj; ++jj)
4532 for (int ii=-ri; ii<=ri; ++ii)
4533 if (mask[rj+jj][ri+ii] && res_sup(i+ii, j+jj)<v)
4534 res_sup(i+ii,j+jj)=0.0f;
4535 }
4536
4537 if (!output_response_mask&&!signed_response)
4538 return res_sup;
4539 unsigned np = 2;
4540 if (output_response_mask&&signed_response)
4541 np = 3;
4542 //response plane and mask plane and/or signed response
4543
4544 vil_image_view<float> res_mask(ni, nj, np);
4545 res_mask.fill(0.0f);
4546 // always copy response to plane 0
4547 for (int j = rj; j<(nj-rj); j++)
4548 for (int i = ri; i<(ni-ri); i++)
4549 {
4550 float rv = res_sup(i,j);
4551 res_mask(i,j,0)=rv;
4552 }
4553 // copy mask plane to plane 1 (or not)
4554 if (output_response_mask) {
4555 for (int j = rj; j<(nj-rj); j++)
4556 for (int i = ri; i<(ni-ri); i++)
4557 {
4558 float rv = res_sup(i,j);
4559 if (!rv)
4560 continue;
4561 for (int jj=-rj; jj<=rj; ++jj)
4562 for (int ii=-ri; ii<=ri; ++ii)
4563 if (mask[jj+rj][ii+ri])
4564 if (rv>res_mask(i+ii,j+jj,1))
4565 res_mask(i+ii,j+jj,1) = rv;
4566 }
4567 }
4568 // copy signed response to plane 1 if no mask plane, otherwise plane 2
4569 if (signed_response) {
4570 unsigned p = np-1;
4571 for (int j = rj; j<(nj-rj); j++)
4572 for (int i = ri; i<(ni-ri); i++)
4573 res_mask(i,j,p) = temp2(i,j);
4574 }
4575 return res_mask;
4576 }
4577
4578 vil_image_view<float> brip_vil_float_ops::
fast_extrema_rotational(vil_image_view<float> const & input,float lambda0,float lambda1,float theta_interval,bool bright,bool mag_only,bool signed_response,bool scale_invariant,bool non_max_suppress,float cutoff,float theta_init,float theta_end)4579 fast_extrema_rotational(vil_image_view<float> const& input,
4580 float lambda0, float lambda1,
4581 float theta_interval,
4582 bool bright, bool mag_only,
4583 bool signed_response, bool scale_invariant,
4584 bool non_max_suppress,
4585 float cutoff, float theta_init, float theta_end) {
4586 unsigned ni = input.ni(), nj = input.nj();
4587 vil_image_view<float> resp(ni, nj, 3);
4588 resp.fill(0.0f);
4589
4590 // elliptical operator has 180 degree rotational symmetry, so only the angles in the range [0,180] matter
4591 float theta = theta_init; unsigned theta_i = 1;
4592 for (; theta < theta_end; theta += theta_interval, theta_i++)
4593 {
4594 vil_image_view<float> temp =
4595 fast_extrema(input, lambda0, lambda1, theta, bright,
4596 mag_only, false, signed_response,
4597 scale_invariant, non_max_suppress, cutoff);//scale_invariant, false, cutoff);
4598 for (unsigned j = 0; j<nj; ++j)
4599 for (unsigned i = 0; i<ni; ++i) {
4600 float v = temp(i,j);
4601 if (std::fabs(v)>std::fabs(resp(i,j,0)))
4602 resp(i,j,0) = v;
4603 }
4604 }
4605 return resp;
4606 }
4607