contrib/brl/bseg/brip/brip_vil_float_ops.cxx
Go to the documentation of this file.
00001 #include "brip_vil_float_ops.h"
00002 //:
00003 // \file
00004 
00005 #include <vcl_cassert.h>
00006 #include <vcl_fstream.h>
00007 #include <vcl_complex.h>
00008 #include <vcl_limits.h>
00009 #include <vgl/vgl_box_2d.h>
00010 #include <vgl/vgl_point_2d.h>
00011 #include <vgl/vgl_vector_2d.h>
00012 #include <vgl/vgl_polygon_scan_iterator.h>
00013 #include <vul/vul_timer.h>
00014 #include <vbl/vbl_array_1d.h>
00015 #include <vbl/vbl_array_2d.h>
00016 #include <vbl/vbl_bounding_box.h>
00017 #include <vnl/vnl_numeric_traits.h>
00018 #include <vnl/vnl_math.h>
00019 #include <vnl/vnl_double_2x3.h>
00020 #include <vnl/algo/vnl_fft_prime_factors.h>
00021 #include <vnl/algo/vnl_svd.h>
00022 
00023 #include <vil/vil_pixel_format.h>
00024 #include <vil/vil_transpose.h>
00025 #include <vil/vil_convert.h>
00026 #include <vil/vil_save.h>
00027 #include <vil/vil_new.h>
00028 #include <vil/vil_math.h>
00029 #include <vil/algo/vil_convolve_1d.h>
00030 
00031 #include <vsol/vsol_box_2d.h>
00032 #include <vsol/vsol_polygon_2d_sptr.h>
00033 #include <vsol/vsol_polygon_2d.h>
00034 #include <bsol/bsol_algs.h>
00035 #include <bsta/bsta_histogram.h>
00036 #include <bsta/bsta_joint_histogram.h>
00037 #include <brip/brip_roi.h>
00038 
00039 // === Local utility functions ===
00040 
00041 //: compute normalized cross correlation from the intensity moment sums.
00042 static float cross_corr(double area, double si1, double si2,
00043                         double si1i1, double si2i2, double si1i2,
00044                         float intensity_thresh)
00045 {
00046   if (!area)
00047     return 0.f;
00048   // the mean values
00049   double u1 = si1/area, u2 = si2/area;
00050   if (u1<intensity_thresh||u2<intensity_thresh)
00051     return -1.f;
00052   double neu = si1i2 - area*u1*u2;
00053   double sd1 = vcl_sqrt(vcl_fabs(si1i1-area*u1*u1)),
00054     sd2 = vcl_sqrt(vcl_fabs(si2i2-area*u2*u2));
00055   if (!neu)
00056     return 0.f;
00057   if (!sd1||!sd2) {
00058     if (neu>0)
00059       return 1.f;
00060     else
00061       return -1.f;
00062   }
00063   double den = sd1*sd2;
00064   return float(neu/den);
00065 }
00066 
00067 //------------------------------------------------------------
00068 //:  Convolve with a kernel
00069 //   It's assumed that the kernel is square with odd dimensions
00070 vil_image_view<float>
00071 brip_vil_float_ops::convolve(vil_image_view<float> const& input,
00072                              vbl_array_2d<float> const& kernel)
00073 {
00074   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
00075   int kw = kernel.cols(); // kh = kernel.rows();
00076   // add a check for kernels that are not equal dimensions of odd size JLM
00077   int n = (kw-1)/2;
00078   vil_image_view<float> output;
00079   output.set_size(w,h);
00080   for (int y = n; y<(h-n); y++)
00081     for (int x = n; x<(w-n); x++)
00082     {
00083       float accum = 0;
00084       for (int j = -n; j<=n; j++)
00085         for (int i = -n; i<=n; i++)
00086         {
00087           float x1 = input(x+i,y+j);
00088           float x2 = kernel[i+n][j+n];
00089           accum += x1*x2;
00090         }
00091       output(x,y)=accum;
00092     }
00093   brip_vil_float_ops::fill_x_border(output, n, 0.0f);
00094   brip_vil_float_ops::fill_y_border(output, n, 0.0f);
00095   return output;
00096 }
00097 
00098 static void fill_1d_array(vil_image_view<float> const& input,
00099                           int y, float* output)
00100 {
00101   unsigned w = input.ni();
00102   for (unsigned x = 0; x<w; x++)
00103     output[x] = input(x,y);
00104 }
00105 
00106 //: Downsamples the 1-d array by 2 using the Burt-Adelson reduction algorithm.
00107 void brip_vil_float_ops::half_resolution_1d(const float* input, unsigned width,
00108                                             float k0, float k1, float k2,
00109                                             float* output)
00110 {
00111   float w[5];
00112   int n = 0;
00113   for (; n<5; n++)
00114     w[n]=input[n];
00115   output[0]=k0*w[0]+ 2.0f*(k1*w[1] + k2*w[2]); // reflect at boundary
00116   for (unsigned x = 1; x<width; ++x)
00117   {
00118     output[x]=k0*w[2]+ k1*(w[1]+w[3]) + k2*(w[0]+w[4]);
00119     //shift the window, w, over by two pixels
00120     w[0] = w[2];       w[1] = w[3];     w[2] = w[4];
00121     //handle the boundary conditions
00122     if (x+2<width)
00123       w[3] = input[n++], w[4] = input[n++];
00124     else
00125       w[3] = w[1], w[4] = w[0];
00126   }
00127 }
00128 
00129 //: Downsamples the image by 2 using the Burt-Adelson reduction algorithm.
00130 // Convolution with a 5-point kernel [(0.5-ka)/2, 0.25, ka, 0.25, (0.5-ka)/2]
00131 // ka = 0.6  maximum decorrelation, wavelet for image compression.
00132 // ka = 0.5  linear interpolation,
00133 // ka = 0.4  Gaussian filter
00134 // ka = 0.359375 min aliasing, wider than Gaussian
00135 // The image sizes are related by: output_dimension = (input_dimension +1)/2.
00136 vil_image_view<float>
00137 brip_vil_float_ops::half_resolution(vil_image_view<float> const& input,
00138                                     float filter_coef)
00139 {
00140   vul_timer t;
00141   float k0 = filter_coef, k1 = 0.25f*filter_coef, k2 = 0.5f*(0.5f-filter_coef);
00142   unsigned w = input.ni(), h = input.nj();
00143   int half_w =(w+1)/2, half_h = (h+1)/2;
00144   vil_image_view<float> output;
00145   output.set_size(half_w, half_h);
00146   // Generate input/output arrays
00147   int n = 0;
00148   float* in0 = new float[w];  float* in1 = new float[w];
00149   float* in2 = new float[w];  float* in3 = new float[w];
00150   float* in4 = new float[w];
00151 
00152   float* out0 = new float[half_w];  float* out1 = new float[half_w];
00153   float* out2 = new float[half_w];  float* out3 = new float[half_w];
00154   float* out4 = new float[half_w];
00155   // Initialize arrays
00156   fill_1d_array(input, n++, in0);   fill_1d_array(input, n++, in1);
00157   fill_1d_array(input, n++, in2);   fill_1d_array(input, n++, in3);
00158   fill_1d_array(input, n++, in4);
00159 
00160   // downsample initial arrays
00161   brip_vil_float_ops::half_resolution_1d(in0, half_w, k0, k1, k2, out0);
00162   brip_vil_float_ops::half_resolution_1d(in1, half_w, k0, k1, k2, out1);
00163   brip_vil_float_ops::half_resolution_1d(in2, half_w, k0, k1, k2, out2);
00164   brip_vil_float_ops::half_resolution_1d(in3, half_w, k0, k1, k2, out3);
00165   brip_vil_float_ops::half_resolution_1d(in4, half_w, k0, k1, k2, out4);
00166   int x=0, y;
00167   // do the first output line
00168   for (;x<half_w;x++)
00169     output(x,0)= k0*out0[x]+ 2.0f*(k1*out1[x]+k2*out2[x]);
00170   // normal lines
00171   for (y=1; y<half_h; y++)
00172   {
00173     for (x=0; x<half_w; x++)
00174       output(x,y) = k0*out2[x]+ k1*(out1[x]+out3[x]) + k2*(out0[x]+out4[x]);
00175     //shift the neighborhood down two lines
00176     float* temp0 = out0;
00177     float* temp1 = out1;
00178     out0 = out2;  out1 = out3;  out2 = out4;
00179     out3 = temp0; out4 = temp1;//reflect values
00180     //test border condition
00181     if (y<half_h-2)
00182     {
00183       //normal processing, so don't reflect
00184       fill_1d_array(input, n++, in3);
00185       fill_1d_array(input, n++, in4);
00186       brip_vil_float_ops::half_resolution_1d(in3, half_w, k0, k1, k2, out3);
00187       brip_vil_float_ops::half_resolution_1d(in4, half_w, k0, k1, k2, out4);
00188     }
00189   }
00190   delete [] in0;  delete [] in1; delete [] in2;
00191   delete [] in3;  delete [] in4;
00192   delete [] out0;  delete [] out1; delete [] out2;
00193   delete [] out3;  delete [] out4;
00194 #ifdef DEBUG
00195   vcl_cout << "\nDownsample a "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
00196 #endif
00197   return output;
00198 }
00199 
00200 void brip_vil_float_ops::double_resolution_1d(const float* input, const unsigned n_input,
00201                                               const float k0, const float k1,
00202                                               const float k2, float* output)
00203 {
00204   float w[3];
00205   unsigned i = 0;
00206   w[1]=input[i]; w[2]=input[i++];
00207   w[0]=w[2];
00208   for (unsigned c = 0; c<2*n_input; c+=2)
00209   {
00210     output[c] = k0*w[1] + k2*(w[0]+w[2]);
00211     output[c+1] = k1*(w[1]+w[2]);
00212     w[0]=w[1];
00213     w[1]=w[2];
00214     if (c<2*(n_input-2))
00215       w[2]=input[i++];
00216     else
00217       w[2]=w[0];
00218   }
00219 }
00220 
00221 //: interpolates the input using the Bert-Adelson algorithm
00222 vil_image_view<float>
00223 brip_vil_float_ops::double_resolution(vil_image_view<float> const& input,
00224                                       float filter_coef)
00225 {
00226   unsigned ni_in = input.ni();
00227   unsigned nj_in = input.nj();
00228   unsigned ni_out = 2*ni_in;
00229   unsigned nj_out = 2*nj_in;
00230   vil_image_view<float> out(ni_out, nj_out);
00231   float* input_1d = new float[ni_in];
00232 
00233   // An interpolation neighborhood of three lines
00234   float* output0 = new float[ni_out];
00235   float* output1 = new float[ni_out];
00236   float* output2 = new float[ni_out];
00237 
00238   // The filter coefficients
00239   float k0 = filter_coef*2.0f;
00240   float k1 = 0.5f;
00241   float k2 = 0.5f-filter_coef;
00242 
00243   // initialize
00244   unsigned i = 0;
00245   fill_1d_array(input, i++, input_1d);
00246   brip_vil_float_ops::double_resolution_1d(input_1d, ni_in, k0, k1, k2, output1);
00247   fill_1d_array(input, i++, input_1d);
00248   brip_vil_float_ops::double_resolution_1d(input_1d, ni_in, k0, k1, k2, output2);
00249   for (unsigned k = 0; k<ni_out; ++k)
00250     output0[k]=output2[k];
00251   for (unsigned r = 0; r<nj_out; r+=2)
00252   {
00253     unsigned rp = r+1;
00254     for (unsigned c=0; c<ni_out; ++c)
00255     {
00256       out(c, r) = k0*output1[c] + k2*(output0[c]+output2[c]);
00257       out(c, rp) = k1*(output1[c]+output2[c]);
00258     }
00259     float* next = output0;
00260     output0 = output1;
00261     output1 = output2;
00262     output2 = next;
00263     if (r<nj_out-4)
00264     {
00265       fill_1d_array(input, i++, input_1d);
00266       brip_vil_float_ops::double_resolution_1d(input_1d, ni_in,
00267                                                k0, k1, k2, output2);
00268     }
00269     else
00270       for (unsigned k = 0; k<ni_out; ++k)
00271         output2[k]=output0[k];
00272   }
00273   delete [] input_1d;
00274   delete [] output0;
00275   delete [] output1;
00276   delete [] output2;
00277   return out;
00278 }
00279 
00280 static double brip_vil_gaussian(double x, double sigma)
00281 {
00282   double x_on_sigma = x / sigma;
00283   return (double)vcl_exp(- x_on_sigma * x_on_sigma / 2);
00284 }
00285 
00286 unsigned brip_vil_float_ops::gaussian_radius(const double sigma,
00287                                              const double fuzz)
00288 {
00289   unsigned radius = 0;
00290   while ((float)brip_vil_gaussian((double)radius, sigma) > fuzz) ++radius;
00291   return radius;
00292 }
00293 
00294 //: generate a 1-d Gaussian kernel  fuzz=0.02 is a good value
00295 static void brip_1d_gaussian_kernel(double sigma, double fuzz,
00296                                     unsigned& radius, double*& kernel)
00297 {
00298   radius = brip_vil_float_ops::gaussian_radius(sigma, fuzz);
00299 
00300   kernel = new double[2*radius + 1];
00301   if (!radius)
00302   {
00303     kernel[0]=1;
00304     return;
00305   }
00306   for (unsigned i=0; i<=radius; ++i)
00307     kernel[radius+i] = kernel[radius-i] = brip_vil_gaussian(double(i), sigma);
00308   double sum = 0;
00309   for (unsigned i= 0; i <= 2*radius; ++i)
00310     sum += kernel[i];                           // find integral of weights
00311   for (unsigned i= 0; i <= 2*radius; ++i)
00312     kernel[i] /= sum;                           // normalize by integral
00313 }
00314 
00315 vil_image_view<float>
00316 brip_vil_float_ops::gaussian(vil_image_view<float> const& input, float sigma,
00317                              vcl_string const& boundary_condition,
00318                              float fill)
00319 {
00320   vil_convolve_boundary_option option=vil_convolve_ignore_edge;
00321   if(boundary_condition == "zeros")
00322     option = vil_convolve_zero_extend;
00323   else if(boundary_condition == "const")
00324     option = vil_convolve_constant_extend;
00325   else if(boundary_condition == "periodic")
00326     option = vil_convolve_periodic_extend;
00327   else if(boundary_condition == "reflect")
00328     option = vil_convolve_reflect_extend;
00329 
00330   unsigned ni = input.ni(), nj = input.nj();
00331   unsigned np = input.nplanes();
00332   vil_image_view<float> dest(ni, nj, np);
00333   dest.fill(fill);
00334   unsigned r;
00335   double* ker;
00336   brip_1d_gaussian_kernel(sigma, 0.02, r, ker);
00337   for (unsigned p = 0; p<np; ++p) {
00338     vil_image_view<float> input_temp(ni, nj);
00339     vil_image_view<float> out_temp(ni, nj);
00340     vil_image_view<float> work(ni, nj);
00341     work.fill(fill);
00342     for (unsigned j = 0; j<nj; ++j)
00343       for (unsigned i = 0; i<ni; ++i)
00344         input_temp(i,j) = input(i,j,p);
00345 
00346     // filter horizontal
00347     int ksize = 2*r + 1 ;
00348     float accum=0.0f;
00349     vil_convolve_1d(input_temp, work, ker + ksize/2,
00350                     -ksize/2, r, accum,
00351                     option, option);
00352     // filter vertical
00353     vil_image_view<float> work_t = vil_transpose(work);
00354     vil_image_view<float> out_temp_t = vil_transpose(dest);
00355     vil_convolve_1d(work_t, out_temp_t, ker+ ksize/2,
00356                     -ksize/2, r, accum,
00357                     option, option);
00358 
00359     vil_image_view<float> plane = vil_transpose(out_temp_t);
00360     for (unsigned j = 0; j<nj; ++j)
00361       for (unsigned i = 0; i<ni; ++i)
00362         dest(i,j,p) = plane(i,j);
00363   }
00364   delete ker;
00365   return dest;
00366 }
00367 
00368 vil_image_view<float>
00369 brip_vil_float_ops::absolute_value(vil_image_view<float> const& input)
00370 {
00371   unsigned ni = input.ni(), nj = input.nj();
00372   unsigned np = input.nplanes();
00373   vil_image_view<float> dest(ni, nj, np);
00374   for (unsigned p = 0; p<np; ++p)
00375     for (unsigned j = 0; j<nj; ++j)
00376       for (unsigned i = 0; i<ni; ++i)
00377         dest(i,j,p) = vcl_fabs(input(i,j,p));
00378   return dest;
00379 }
00380 
00381 #ifdef VIL_CONVOLVE_WITH_MASK_EXISTS // TODO
00382 vil_image_view<float>
00383 brip_vil_float_ops::gaussian(vil_image_view<float> const& input,
00384                              float sigma,
00385                              vil_image_view<float> const& mask)
00386 {
00387   vil_image_view<float> dest(input.ni(), input.nj());
00388 
00389   int r;
00390   double* ker;
00391   brip_1d_gaussian_kernel(sigma, 0.02, r, ker);
00392   vil_image_view<float> work(input.ni(), input.nj());
00393   work.deep_copy(input);
00394   // filter horizontal
00395   int ksize = 2*r + 1 ;
00396   float accum=0.0f;
00397   vil_convolve_1d(input, work, mask, ker + ksize/2,
00398                   -ksize/2, r, accum,
00399                   vil_convolve_trim,
00400                   vil_convolve_trim);
00401 
00402   // filter vertical
00403   vil_image_view<float> work_t = vil_transpose(work);
00404   vil_image_view<float> dest_t = vil_transpose(dest);
00405   vil_convolve_1d(work_t, dest_t, vil_transpose(mask), ker+ ksize/2,
00406                   -ksize/2, r, accum,
00407                   vil_convolve_constant_extend,
00408                   vil_convolve_constant_extend);
00409 
00410   delete ker;
00411   return dest;
00412 }
00413 
00414 #endif // VIL_CONVOLVE_WITH_MASK_EXISTS
00415 
00416 //-------------------------------------------------------------------
00417 //: Determine if the center of a (2n+1)x(2n+1) neighborhood is a local maximum
00418 //
00419 bool brip_vil_float_ops::
00420 local_maximum(vbl_array_2d<float> const& neighborhood, int n, float& value)
00421 {
00422   bool local_max = true;
00423   value = 0;
00424   float center = neighborhood[n][n];
00425   for (int y = -n; y<=n; y++)
00426     for (int x = -n; x<=n; x++)
00427       local_max = local_max&&(neighborhood[y+n][x+n]<=center);
00428   if (!local_max)
00429     return false;
00430   value = center;
00431   return true;
00432 }
00433 
00434 //-------------------------------------------------------------------
00435 // Interpolate the sub-pixel position of a neighborhood using a
00436 // second order expansion on a 3x3 sub-neighborhood. Return the
00437 // offset to the maximum, i.e. x=x0+dx, y = y0+dy. The design is
00438 // similar to the droid counterpoint by fsm, which uses the Beaudet Hessian
00439 //
00440 void brip_vil_float_ops::
00441 interpolate_center(vbl_array_2d<float>const& neighborhood, float& dx, float& dy)
00442 {
00443   dx = 0; dy=0;
00444   // extract the neighborhood
00445   float n_m1_m1 = neighborhood[0][0];
00446   float n_m1_0 = neighborhood[0][1];
00447   float n_m1_1 = neighborhood[0][2];
00448   float n_0_m1 = neighborhood[1][0];
00449   float n_0_0 = neighborhood[1][1];
00450   float n_0_1 = neighborhood[1][2];
00451   float n_1_m1 = neighborhood[2][0];
00452   float n_1_0 = neighborhood[2][1];
00453   float n_1_1 = neighborhood[2][2];
00454 
00455   // Compute the 2nd order quadratic coefficients
00456   //      1/6 * [ -1  0 +1 ]
00457   // Ix =       [ -1  0 +1 ]
00458   //            [ -1  0 +1 ]
00459   float Ix =(-n_m1_m1+n_m1_1-n_0_m1+n_0_1-n_1_m1+n_1_1)/6.0f;
00460   //      1/6 * [ -1 -1 -1 ]
00461   // Iy =       [  0  0  0 ]
00462   //            [ +1 +1 +1 ]
00463   float Iy =(-n_m1_m1-n_m1_0-n_m1_1+n_1_m1+n_1_0+n_1_1)/6.0f;
00464   //      1/3 * [ +1 -2 +1 ]
00465   // Ixx =      [ +1 -2 +1 ]
00466   //            [ +1 -2 +1 ]
00467   float Ixx = ((n_m1_m1+n_0_m1+n_1_m1+n_m1_1+n_0_1+n_1_1)
00468                -2.0f*(n_m1_0+n_0_0+n_1_0))/3.0f;
00469   //      1/4 * [ +1  0 -1 ]
00470   // Ixy =      [  0  0  0 ]
00471   //            [ -1  0 +1 ]
00472   float Ixy = (n_m1_m1-n_m1_1+n_1_m1+n_1_1)/4.0f;
00473   //      1/3 * [ +1 +1 +1 ]
00474   // Iyy =      [ -2 -2 -2 ]
00475   //            [ +1 +1 +1 ]
00476   float Iyy = ((n_m1_m1+n_m1_0+n_m1_1+n_1_m1+n_1_0+n_1_1)
00477                -2.0f*(n_0_m1 + n_0_0 + n_1_0))/3.0f;
00478   //
00479   // The next bit is to find the extremum of the fitted surface by setting its
00480   // partial derivatives to zero. We need to solve the following linear system :
00481   // Given the fitted surface is
00482   // I(x,y) = Io + Ix x + Iy y + 1/2 Ixx x^2 + Ixy x y + 1/2 Iyy y^2
00483   // we solve for the maximum (x,y),
00484   //
00485   //  [ Ixx Ixy ] [ dx ] + [ Ix ] = [ 0 ]      (dI/dx = 0)
00486   //  [ Ixy Iyy ] [ dy ]   [ Iy ]   [ 0 ]      (dI/dy = 0)
00487   //
00488   float det = Ixx*Iyy - Ixy*Ixy;
00489   // det>0 corresponds to a true local extremum otherwise a saddle point
00490   if (det>0)
00491   {
00492     dx = (Iy*Ixy - Ix*Iyy) / det;
00493     dy = (Ix*Ixy - Iy*Ixx) / det;
00494     // more than one pixel away
00495     if (vcl_fabs(dx) > 1.0 || vcl_fabs(dy) > 1.0)
00496       dx = 0; dy = 0;
00497   }
00498 }
00499 
00500 //---------------------------------------------------------------
00501 // Compute the local maxima of the input on a (2n+1)x(2n+2)
00502 // neighborhood above the given threshold. At each local maximum,
00503 // compute the sub-pixel location, (x_pos, y_pos).
00504 void brip_vil_float_ops::
00505 non_maximum_suppression(vil_image_view<float> const& input,
00506                         int n, float thresh,
00507                         vcl_vector<float>& x_pos,
00508                         vcl_vector<float>& y_pos,
00509                         vcl_vector<float>& value)
00510 {
00511   vul_timer t;
00512   int N = 2*n+1;
00513   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
00514   x_pos.clear();  x_pos.clear();   value.clear();
00515   vbl_array_2d<float> neighborhood(N,N);
00516   for (int y =n; y<h-n; y++)
00517     for (int x = n; x<w-n; x++)
00518     {
00519       //If the center is not above threshold then there is
00520       //no hope
00521       if (input(x,y)<thresh)
00522         continue;
00523       //Fill the neighborhood
00524       for (int i = -n; i<=n; i++)
00525         for (int j = -n; j<=n; j++)
00526           neighborhood.put(j+n,i+n,input(x+i, y+j));
00527       //Check if the center is a local maximum
00528       float dx, dy, max_v;
00529       if (brip_vil_float_ops::local_maximum(neighborhood, n, max_v))
00530       {
00531         //if so sub-pixel interpolate (3x3) and output results
00532         brip_vil_float_ops::interpolate_center(neighborhood, dx, dy);
00533         x_pos.push_back((float)x+dx);
00534         y_pos.push_back((float)y+dy);
00535         value.push_back(max_v);
00536       }
00537     }
00538 #ifdef DEBUG
00539   vcl_cout << "\nCompute non-maximum suppression on a "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
00540 #endif
00541 }
00542 
00543 // -----------------------------------------------------------------
00544 //: Subtract image_1 from image_2.
00545 // Will not operate unless the two input images are the same dimensions
00546 //
00547 vil_image_view<float>
00548 brip_vil_float_ops::difference(vil_image_view<float> const& image_1,
00549                                vil_image_view<float> const& image_2)
00550 {
00551   unsigned w1 = image_1.ni(), h1 = image_1.nj();
00552   unsigned w2 = image_2.ni(), h2 = image_2.nj();
00553   vil_image_view<float> temp(w1, h1);
00554   if (w1!=w2||h1!=h2)
00555   {
00556     vcl_cout << "In brip_vil_float_ops::difference(..) - images are not the same dimensions\n";
00557     return temp;
00558   }
00559   vil_image_view<float> out;
00560   out.set_size(w1, h1);
00561   for (unsigned y = 0; y<h1; y++)
00562     for (unsigned x = 0; x<w1; x++)
00563       out(x,y) = image_2(x,y)-image_1(x,y);
00564   return out;
00565 }
00566 
00567 //: negate an image returning the same pixel type (only greyscale)
00568 vil_image_resource_sptr brip_vil_float_ops::negate(vil_image_resource_sptr const& imgr)
00569 {
00570   vil_image_resource_sptr outr;
00571   if (!imgr)
00572     return outr;
00573 
00574   vil_pixel_format fmt = imgr->pixel_format();
00575   switch (fmt)
00576     {
00577 #define NEGATE_CASE(FORMAT, T) \
00578    case FORMAT: { \
00579     vil_image_view<T> view = imgr->get_copy_view(); \
00580     T mxv = vcl_numeric_limits<T>::max(); \
00581     vil_math_scale_and_offset_values(view, -1.0, mxv); \
00582     outr = vil_new_image_resource_of_view(view);  \
00583     break; \
00584                 }
00585       NEGATE_CASE(VIL_PIXEL_FORMAT_BYTE, vxl_byte);
00586       NEGATE_CASE(VIL_PIXEL_FORMAT_UINT_32, vxl_uint_32);
00587       NEGATE_CASE(VIL_PIXEL_FORMAT_UINT_16, vxl_uint_16);
00588       NEGATE_CASE(VIL_PIXEL_FORMAT_INT_16, vxl_int_16);
00589       NEGATE_CASE(VIL_PIXEL_FORMAT_FLOAT, float);
00590       NEGATE_CASE(VIL_PIXEL_FORMAT_DOUBLE, double);
00591 #undef NEGATE_CASE
00592     default:
00593       vcl_cout << "Unknown image format\n";
00594     }
00595   return outr;
00596 }
00597 
00598 vil_image_view<float>
00599 brip_vil_float_ops::threshold(vil_image_view<float> const & image,
00600                               const float thresh, const float level)
00601 {
00602   vil_image_view<float> out;
00603   unsigned w = image.ni(), h = image.nj();
00604   out.set_size(w, h);
00605   for (unsigned y = 0; y<h; y++)
00606     for (unsigned x = 0; x<w; x++)
00607     {
00608       if (image(x,y)>thresh)
00609         out(x,y) = level;
00610       else
00611         out(x,y) = 0;
00612     }
00613   return out;
00614 }
00615 
00616 vil_image_view<float>
00617 brip_vil_float_ops::abs_clip_to_level(vil_image_view<float> const& image,
00618                                       float thresh, float level)
00619 {
00620   vil_image_view<float> out;
00621   unsigned w = image.ni(), h = image.nj();
00622   out.set_size(w, h);
00623   for (unsigned y = 0; y<h; y++)
00624     for (unsigned x = 0; x<w; x++)
00625     {
00626       if (vcl_fabs(image(x,y))>thresh)
00627         out(x,y) = level;
00628       else
00629         out(x,y) = image(x,y);
00630     }
00631   return out;
00632 }
00633 
00634 vil_image_view<float> brip_vil_float_ops::average_NxN(vil_image_view<float> const & img, int N)
00635 {
00636   vil_image_view<float> result;
00637   unsigned w = img.ni(), h = img.nj();
00638   result.set_size (w, h);
00639 
00640   vbl_array_2d <float> averaging_filt(N, N);
00641   averaging_filt.fill( float(1.00/double(N*N)) );
00642   result = brip_vil_float_ops::convolve(img, averaging_filt);
00643   return result;
00644 }
00645 
00646 //----------------------------------------------------------------
00647 //: Compute the gradient of the input, use a 3x3 mask
00648 // \verbatim
00649 //         1  |-1  0  1|         1  |-1 -1 -1|
00650 //   Ix = --- |-1  0  1|   Iy = --- | 0  0  0|
00651 //         6  |-1  0  1|         6  | 1  1  1|
00652 // \endverbatim
00653 // Larger masks are computed by pre-convolving with a Gaussian
00654 //
00655 void brip_vil_float_ops::gradient_3x3(vil_image_view<float> const& input,
00656                                       vil_image_view<float>& grad_x,
00657                                       vil_image_view<float>& grad_y)
00658 {
00659   vul_timer t;
00660   unsigned w = input.ni(), h = input.nj();
00661   float scale = 1.0f/6.0f;
00662   for (unsigned y = 1; y+1<h; ++y)
00663     for (unsigned x = 1; x+1<w; ++x)
00664     {
00665       float gx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
00666                 -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
00667       float gy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
00668                 -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
00669       grad_x(x,y) = scale*gx;
00670       grad_y(x,y) = scale*gy;
00671     }
00672   brip_vil_float_ops::fill_x_border(grad_x, 1, 0.0f);
00673   brip_vil_float_ops::fill_y_border(grad_x, 1, 0.0f);
00674   brip_vil_float_ops::fill_x_border(grad_y, 1, 0.0f);
00675   brip_vil_float_ops::fill_y_border(grad_y, 1, 0.0f);
00676 #ifdef DEBUG
00677   vcl_cout << "\nCompute Gradient in " << t.real() << " msecs.\n";
00678 #endif
00679 }
00680 
00681 void brip_vil_float_ops::gradient_mag_3x3(vil_image_view<float> const& input,
00682                                           vil_image_view<float>& mag)
00683 {
00684   unsigned w = input.ni(), h = input.nj();
00685   float scale = 1.0f/6.0f;
00686   for (unsigned y = 1; y+1<h; ++y)
00687     for (unsigned x = 1; x+1<w; ++x)
00688     {
00689       float gx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
00690                 -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
00691       float gy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
00692                 -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
00693       mag(x,y) = scale*vcl_sqrt(gx*gx+gy*gy);
00694     }
00695   brip_vil_float_ops::fill_x_border(mag, 1, 0.0f);
00696   brip_vil_float_ops::fill_y_border(mag, 1, 0.0f);
00697 }
00698 
00699 void brip_vil_float_ops::
00700 gradient_mag_comp_3x3(vil_image_view<float> const& input,
00701                       vil_image_view<float>& mag,
00702                       vil_image_view<float>& gx,
00703                       vil_image_view<float>& gy)
00704 {
00705   unsigned w = input.ni(), h = input.nj();
00706   float scale = 1.0f/6.0f;
00707   for (unsigned y = 1; y+1<h; ++y)
00708     for (unsigned x = 1; x+1<w; ++x)
00709     {
00710       float ggx = input(x+1,y-1)+input(x+1,y)+input(x+1,y-1)
00711                  -input(x-1,y-1)-input(x-1,y)-input(x-1,y-1);
00712       float ggy = input(x+1,y+1)+input(x,y+1)+input(x-1,y+1)
00713                  -input(x+1,y-1)-input(x,y-1)-input(x-1,y-1);
00714       mag(x,y) = scale*vcl_sqrt(ggx*ggx+ggy*ggy);
00715       gx(x,y) = ggx*scale;
00716       gy(x,y) = ggy*scale;
00717     }
00718   brip_vil_float_ops::fill_x_border(mag, 1, 0.0f);
00719   brip_vil_float_ops::fill_y_border(mag, 1, 0.0f);
00720   brip_vil_float_ops::fill_x_border(gx, 1, 0.0f);
00721   brip_vil_float_ops::fill_y_border(gy, 1, 0.0f);
00722 }
00723 
00724 //----------------------------------------------------------------
00725 //: Compute the Hessian of the input, use a 3x3 mask
00726 // \verbatim
00727 //          1 | 1  -2  1|          1 |  1  1  1|         1  | 1  0 -1|
00728 //   Ixx = ---| 1  -2  1|   Iyy = ---| -2 -2 -2|  Ixy = --- | 0  0  0|
00729 //          3 | 1  -2  1|          3 |  1  1  1|         4  |-1  0  1|
00730 // \endverbatim
00731 // Larger masks are computed by pre-convolving with a Gaussian
00732 //
00733 void brip_vil_float_ops::hessian_3x3(vil_image_view<float> const& input,
00734                                      vil_image_view<float>& Ixx,
00735                                      vil_image_view<float>& Ixy,
00736                                      vil_image_view<float>& Iyy)
00737 {
00738   vul_timer t;
00739   unsigned w = input.ni(), h = input.nj();
00740   for (unsigned y = 1; y+1<h; ++y)
00741     for (unsigned x = 1; x+1<w; ++x)
00742     {
00743       float xx = input(x-1,y-1)+input(x-1,y)+input(x+1,y+1)+
00744                  input(x+1,y-1)+input(x+1,y)+input(x+1,y+1)-
00745                  2.0f*(input(x,y-1)+input(x,y)+input(x,y+1));
00746 
00747       float xy = (input(x-1,y-1)+input(x+1,y+1))-
00748                  (input(x-1,y+1)+input(x+1,y-1));
00749 
00750       float yy = input(x-1,y-1)+input(x,y-1)+input(x+1,y-1)+
00751                  input(x-1,y+1)+input(x,y+1)+input(x+1,y+1)-
00752                  2.0f*(input(x-1,y)+input(x,y)+input(x+1,y));
00753 
00754       Ixx(x,y) = xx/3.0f;
00755       Ixy(x,y) = xy/4.0f;
00756       Iyy(x,y) = yy/3.0f;
00757     }
00758   brip_vil_float_ops::fill_x_border(Ixx, 1, 0.0f);
00759   brip_vil_float_ops::fill_y_border(Ixx, 1, 0.0f);
00760   brip_vil_float_ops::fill_x_border(Ixy, 1, 0.0f);
00761   brip_vil_float_ops::fill_y_border(Ixy, 1, 0.0f);
00762   brip_vil_float_ops::fill_x_border(Iyy, 1, 0.0f);
00763   brip_vil_float_ops::fill_y_border(Iyy, 1, 0.0f);
00764 #ifdef DEBUG
00765   vcl_cout << "\nCompute a hessian matrix "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
00766 #endif
00767 }
00768 
00769 vil_image_view<float>
00770 brip_vil_float_ops::beaudet(vil_image_view<float> const& Ixx,
00771                             vil_image_view<float> const& Ixy,
00772                             vil_image_view<float> const& Iyy,
00773                             bool determinant)
00774 {
00775   unsigned w = Ixx.ni(), h = Ixx.nj();
00776   vil_image_view<float> output;
00777   output.set_size(w, h);
00778   for (unsigned y = 0; y<h; y++)
00779     for (unsigned x = 0; x<w; x++)
00780     {
00781       float xx = Ixx(x,y), xy = Ixy(x,y), yy = Iyy(x,y);
00782 
00783       //compute eigenvalues for experimentation
00784       float det = xx*yy-xy*xy;
00785       float tr = xx+yy;
00786       if (determinant)
00787         output(x,y) = det;
00788       else
00789         output(x,y) = tr;
00790     }
00791   return output;
00792 }
00793 
00794 //----------------------------------------------------------------
00795 //: $Ix\cdot Ix^t$ gradient matrix elements
00796 // That is,
00797 // \verbatim
00798 //                        _                           _
00799 //                       | (dI/dx)^2    (dI/dx)(dI/dy) |
00800 //                       |                             |
00801 //  A = Sum(neighborhood)|                             |
00802 //                       |(dI/dx)(dI/dy)   (dI/dx)^2   |
00803 //                       |_                           _|
00804 // \endverbatim
00805 // over a 2n+1 x 2n+1 neighborhood
00806 //
00807 void
00808 brip_vil_float_ops::grad_matrix_NxN(vil_image_view<float> const& input,
00809                                     unsigned n,
00810                                     vil_image_view<float>& IxIx,
00811                                     vil_image_view<float>& IxIy,
00812                                     vil_image_view<float>& IyIy)
00813 {
00814   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
00815   int N = (2*n+1)*(2*n+1);
00816   int ni = static_cast<int>(n);
00817   vil_image_view<float> grad_x, grad_y, output;
00818   grad_x.set_size(w,h);
00819   grad_y.set_size(w,h);
00820   output.set_size(w,h);
00821   brip_vil_float_ops::gradient_3x3(input, grad_x, grad_y);
00822   vul_timer t;
00823   for (int y = ni; y<h-ni;y++)
00824     for (int x = ni; x<w-ni;x++)
00825     {
00826       float xx=0, xy=0, yy=0;
00827       for (int i = -ni; i<=ni; i++)
00828         for (int j = -ni; j<=ni; j++)
00829         {
00830           float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
00831           xx += gx*gx;
00832           xy += gx*gy;
00833           yy += gy*gy;
00834         }
00835       IxIx(x,y) = xx/(float)N;
00836       IxIy(x,y) = xy/(float)N;
00837       IyIy(x,y) = yy/(float)N;
00838     }
00839   brip_vil_float_ops::fill_x_border(IxIx, ni, 0.0f);
00840   brip_vil_float_ops::fill_y_border(IxIx, ni, 0.0f);
00841   brip_vil_float_ops::fill_x_border(IxIy, ni, 0.0f);
00842   brip_vil_float_ops::fill_y_border(IxIy, ni, 0.0f);
00843   brip_vil_float_ops::fill_x_border(IyIy, ni, 0.0f);
00844   brip_vil_float_ops::fill_y_border(IyIy, ni, 0.0f);
00845 #ifdef DEBUG
00846   vcl_cout << "\nCompute a gradient matrix "<< w <<" x " << h << " image in "<< t.real() << " msecs.\n";
00847 #endif
00848 }
00849 
00850 vil_image_view<float> brip_vil_float_ops::
00851 trace_grad_matrix_NxN(vil_image_view<float> const& input, unsigned n)
00852 {
00853   unsigned ni = input.ni(), nj = input.nj();
00854   vil_image_view<float> IxIx;
00855   vil_image_view<float> IxIy;
00856   vil_image_view<float> IyIy;
00857   vil_image_view<float> tr;
00858   IxIx.set_size(ni, nj);   IxIy.set_size(ni, nj);   IyIy.set_size(ni, nj);
00859   tr.set_size(ni, nj);
00860   brip_vil_float_ops::grad_matrix_NxN(input, n, IxIx, IxIy, IyIy);
00861   vil_math_image_sum<float, float, float>(IxIx, IyIy, tr);
00862   return tr;
00863 }
00864 
00865 vil_image_view<float>
00866 brip_vil_float_ops::harris(vil_image_view<float> const& IxIx,
00867                            vil_image_view<float> const& IxIy,
00868                            vil_image_view<float> const& IyIy,
00869                            double scale)
00870 {
00871   unsigned w = IxIx.ni(), h = IxIx.nj();
00872   float norm = 1e-3f; // Scale the output to values in the 10->1000 range
00873   vil_image_view<float> output;
00874   output.set_size(w, h);
00875   for (unsigned y = 0; y<h; y++)
00876     for (unsigned x = 0; x<w; x++)
00877     {
00878       float xx = IxIx(x,y), xy = IxIy(x,y), yy = IyIy(x,y);
00879       float det = xx*yy-xy*xy, trace = xx+yy;
00880       output(x,y) = float(det - scale*trace*trace)*norm;
00881     }
00882   return output;
00883 }
00884 
00885 //----------------------------------------------------------------
00886 // Compute the sqrt of the product of the eigenvalues of the
00887 // gradient matrix over a 2n+1 x 2n+1 neighborhood
00888 // That is,
00889 // \verbatim
00890 //                        _                           _
00891 //                       | (dI/dx)^2    (dI/dx)(dI/dy) |
00892 //                       |                             |
00893 //  A = Sum(neighborhood)|                             |
00894 //                       |(dI/dx)(dI/dy)   (dI/dx)^2   |
00895 //                       |_                           _|
00896 // \endverbatim
00897 // The output image is sqrt(lamba_1*lambda_2) where lambda_i are the eigenvalues
00898 //
00899 vil_image_view<float>
00900 brip_vil_float_ops::sqrt_grad_singular_values(vil_image_view<float>& input,
00901                                               int n)
00902 {
00903   int N = (2*n+1)*(2*n+1);
00904   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
00905   vil_image_view<float> grad_x, grad_y, output;
00906   grad_x.set_size(w,h);
00907   grad_y.set_size(w,h);
00908   output.set_size(w,h);
00909   brip_vil_float_ops::gradient_3x3(input, grad_x, grad_y);
00910   vul_timer t;
00911   for (int y = n; y<h-n;y++)
00912     for (int x = n; x<w-n;x++)
00913     {
00914       float IxIx=0, IxIy=0, IyIy=0;
00915       for (int i = -n; i<=n; i++)
00916         for (int j = -n; j<=n; j++)
00917         {
00918           float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
00919           IxIx += gx*gx;
00920           IxIy += gx*gy;
00921           IyIy += gy*gy;
00922         }
00923       float det = (IxIx*IyIy-IxIy*IxIy)/(float)N;
00924       output(x,y)=vcl_sqrt(vcl_fabs(det));
00925     }
00926   brip_vil_float_ops::fill_x_border(output, n, 0.0f);
00927   brip_vil_float_ops::fill_y_border(output, n, 0.0f);
00928 #ifdef DEBUG
00929   vcl_cout << "\nCompute sqrt(sigma0*sigma1) in" << t.real() << " msecs.\n";
00930 #endif
00931   return output;
00932 }
00933 
00934 vil_image_view<float> brip_vil_float_ops::
00935 max_scale_trace(vil_image_view<float> input,
00936                 float min_scale, float max_scale, float scale_inc)
00937 {
00938   unsigned ni = input.ni(), nj = input.nj();
00939   vil_image_view<float> tr_max, sc;
00940   tr_max.set_size(ni, nj);
00941   tr_max.fill(0.0f);
00942   sc.set_size(ni, nj);
00943   sc.fill(min_scale);
00944   for (float s = min_scale; s<=max_scale; s+=scale_inc)
00945   {
00946     vil_image_view<float> smooth = brip_vil_float_ops::gaussian(input, s);
00947     unsigned N = static_cast<unsigned>(2.0f*s);
00948     vil_image_view<float> tr =
00949       brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
00950     for (unsigned r = 0; r<nj; ++r)
00951       for (unsigned c = 0; c<ni; ++c)
00952       {
00953         float trv = s*s*tr(c,r);
00954         if (trv>tr_max(c,r))
00955         {
00956           tr_max(c,r) = trv;
00957           sc(c,r) = s;
00958         }
00959       }
00960   }
00961   return sc;
00962 }
00963 
00964 //: 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
00965 vil_image_view<float> brip_vil_float_ops::
00966 max_scale_trace_value(vil_image_view<float> input,
00967                       float min_scale, float max_scale, float scale_inc)
00968 {
00969   unsigned ni = input.ni(), nj = input.nj();
00970   vil_image_view<float> tr_max, sc;
00971   tr_max.set_size(ni, nj);
00972   tr_max.fill(0.0f);
00973   sc.set_size(ni, nj);
00974   sc.fill(min_scale);
00975   for (float s = min_scale; s<=max_scale; s+=scale_inc)
00976   {
00977     vil_image_view<float> smooth = brip_vil_float_ops::gaussian(input, s);
00978     unsigned N = static_cast<unsigned>(2.0f*s);
00979     vil_image_view<float> tr =
00980       brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
00981     for (unsigned r = 0; r<nj; ++r)
00982       for (unsigned c = 0; c<ni; ++c)
00983       {
00984         float trv = s*s*tr(c,r);
00985         if (trv>tr_max(c,r))
00986         {
00987           tr_max(c,r) = trv;
00988           sc(c,r) = s;
00989         }
00990       }
00991   }
00992   // mask the region where integration region extends outside the image borders
00993   unsigned N = static_cast<unsigned>(5.0f*max_scale);
00994   unsigned njmax = nj-N;
00995   unsigned nimax = ni-N;
00996   for (unsigned r = 0; r<nj; ++r)
00997     for (unsigned c = 0; c<ni; ++c)
00998     {
00999       if (r <= N || r >= njmax || c <= N || c >= nimax)
01000         tr_max(c,r) = 0.0f;
01001     }
01002 
01003   // normalize the trace values
01004   float min_b,max_b;
01005   vil_math_value_range(tr_max,min_b,max_b);
01006   vcl_cout << "in trace max image min value: " << min_b << " max value: " << max_b << vcl_endl;
01007 
01008   vil_image_view<double> tr_normalized, tr_stretched;
01009   vil_convert_stretch_range(tr_max, tr_normalized, 0.0f, 1.0f);
01010   vil_convert_stretch_range(tr_max, tr_stretched, 0.0f, 256.0f);
01011   vil_image_view<float> tr_cast;
01012   vil_convert_cast(tr_stretched, tr_cast);
01013   vil_image_view<vxl_byte> tr_cast_byte;
01014   vil_convert_cast(tr_stretched, tr_cast_byte);
01015   vil_save_image_resource(vil_new_image_resource_of_view(tr_cast_byte), "D:\\projects\\vehicle_rec_on_models\\trace_image.png");
01016 
01017 #if 0 // commented out
01018   // investigate illumination invariance
01019   vil_image_view<float> modified_input = input;
01020   vil_image_view<float> const_image(ni, nj);
01021   const_image.fill(3.0f);
01022   for (unsigned r = 0; r<nj; ++r)
01023     for (unsigned c = 0; c<ni; ++c)
01024     {
01025       if (modified_input(c,r) < 120 || modified_input(c,r) > 95) = modified_input(c,r)*const_image(c,r);
01026     }
01027 
01028   vil_image_view<vxl_byte> mod_cast;
01029   vil_convert_cast(modified_input, mod_cast);
01030   vil_save_image_resource(vil_new_image_resource_of_view(mod_cast), "D:\\projects\\vehicle_rec_on_models\\modified_image.png");
01031 
01032   vil_image_view<float> tr_max2, sc2;
01033   tr_max2.set_size(ni, nj);
01034   tr_max2.fill(0.0f);
01035   sc2.set_size(ni, nj);
01036   sc2.fill(min_scale);
01037   for (float s = min_scale; s<=max_scale; s+=scale_inc)
01038   {
01039     vil_image_view<float> smooth = brip_vil_float_ops::gaussian(modified_input, s);
01040     unsigned N = static_cast<unsigned>(2.0f*s);
01041     vil_image_view<float> tr =
01042       brip_vil_float_ops::trace_grad_matrix_NxN(smooth, N);
01043     for (unsigned r = 0; r<nj; ++r)
01044       for (unsigned c = 0; c<ni; ++c)
01045       {
01046         float trv = s*s*tr(c,r);
01047         if (trv>tr_max2(c,r))
01048         {
01049           tr_max2(c,r) = trv;
01050           sc2(c,r) = s;
01051         }
01052       }
01053   }
01054   // mask the region where integration region extends outside the image borders
01055   for (unsigned r = 0; r<nj; ++r)
01056     for (unsigned c = 0; c<ni; ++c)
01057     {
01058       if (r <= N || r >= njmax || c <= N || c >= nimax)
01059         tr_max2(c,r) = 0.0f;
01060     }
01061 
01062   vil_math_value_range(tr_max2,min_b,max_b);
01063   vcl_cout << "in trace max2 image min value: " << min_b << " max value: " << max_b << vcl_endl;
01064 
01065   // normalize
01066   vil_image_view<double> tr_normalized2;
01067   vil_convert_stretch_range(tr_max2, tr_normalized2, 0.0f, 1.0f);
01068 
01069   vil_image_view<double> difference_image = tr_normalized;
01070   for (unsigned r = 0; r<nj; ++r)
01071     for (unsigned c = 0; c<ni; ++c)
01072     {
01073       difference_image(c,r) = vcl_abs(difference_image(c,r) - tr_normalized2(c,r));
01074     }
01075   double min_bd, max_bd;
01076   vil_math_value_range(difference_image,min_bd,max_bd);
01077   vcl_cout << "in difference image of normalized trace images min value: " << min_bd << " max value: " << max_bd << vcl_endl;
01078   vil_image_view<vxl_byte> dif_cast;
01079   vil_convert_cast(difference_image, dif_cast);
01080   vil_save_image_resource(vil_new_image_resource_of_view(dif_cast), "D:\\projects\\vehicle_rec_on_models\\difference_image.png");
01081 
01082   vil_image_view<double> tr_stretched2;
01083   vil_convert_stretch_range(tr_max2, tr_stretched2, 0.0f, 256.0f);
01084   vil_image_view<float> tr_cast2;
01085   vil_convert_cast(tr_stretched2, tr_cast2);
01086   vil_image_view<vxl_byte> tr_cast2_byte;
01087   vil_convert_cast(tr_stretched2, tr_cast2_byte);
01088   vil_save_image_resource(vil_new_image_resource_of_view(tr_cast2_byte), "D:\\projects\\vehicle_rec_on_models\\trace_image2.png");
01089 
01090   return tr_cast2;
01091 #else // 0
01092   return tr_cast;
01093 #endif // 0
01094 }
01095 
01096 
01097 //---------------------------------------------------------------------
01098 // Lucas-Kanade motion vectors:  Solve for the motion vectors over a
01099 // (2n+1)x(2n+1) neighborhood. The time derivative of intensity is computed
01100 // from the previous_frame. The threshold eliminates small values of
01101 // the product of the time derivative and the motion matrix eigenvalues,
01102 // i.e, |lambda_1*lambda_2*dI/dt|<thresh.  Thus motion is only reported when
01103 // the solution is well-conditioned.
01104 //
01105 void
01106 brip_vil_float_ops::Lucas_KanadeMotion(vil_image_view<float> & current_frame,
01107                                        vil_image_view<float> & previous_frame,
01108                                        int n, double thresh,
01109                                        vil_image_view<float>& vx,
01110                                        vil_image_view<float>& vy)
01111 {
01112   int N = (2*n+1)*(2*n+1);
01113   int w = static_cast<int>(current_frame.ni()), h = static_cast<int>(current_frame.nj());
01114   vil_image_view<float> grad_x, grad_y, diff;
01115   grad_x.set_size(w,h);
01116   grad_y.set_size(w,h);
01117   // compute the gradient vector and the time derivative
01118   brip_vil_float_ops::gradient_3x3(current_frame, grad_x, grad_y);
01119   diff = brip_vil_float_ops::difference(previous_frame, current_frame);
01120   vul_timer t;
01121   // sum the motion terms over the (2n+1)x(2n+1) neighborhood.
01122   for (int y = n; y<h-n;y++)
01123     for (int x = n; x<w-n;x++)
01124     {
01125       float IxIx=0, IxIy=0, IyIy=0, IxIt=0, IyIt=0;
01126       for (int i = -n; i<=n; i++)
01127         for (int j = -n; j<=n; j++)
01128         {
01129           float gx = grad_x(x+i, y+j), gy = grad_y(x+i, y+j);
01130           float dt = diff(x+i, y+j);
01131           IxIx += gx*gx;
01132           IxIy += gx*gy;
01133           IyIy += gy*gy;
01134           IxIt += gx*dt;
01135           IyIt += gy*dt;
01136         }
01137       //Divide by the number of pixels in the neighborhood
01138       IxIx/=N;  IxIy/=N; IyIy/=N; IxIt/=N; IyIt/=N;
01139       float det = float(IxIx*IyIy-IxIy*IxIy);
01140       //Eliminate small motion factors
01141       float dif = diff(x,y);
01142       float motion_factor = vcl_fabs(det*dif);
01143       if (motion_factor<thresh)
01144       {
01145         vx(x,y) = 0.0f;
01146         vy(x,y) = 0.0f;
01147         continue;
01148       }
01149       //solve for the motion vector
01150       vx(x,y) = (IyIy*IxIt-IxIy*IyIt)/det;
01151       vy(x,y) = (-IxIy*IxIt + IxIx*IyIt)/det;
01152     }
01153   brip_vil_float_ops::fill_x_border(vx, n, 0.0f);
01154   brip_vil_float_ops::fill_y_border(vx, n, 0.0f);
01155   brip_vil_float_ops::fill_x_border(vy, n, 0.0f);
01156   brip_vil_float_ops::fill_y_border(vy, n, 0.0f);
01157 #ifdef DEBUG
01158   vcl_cout << "\nCompute Lucas-Kanade in " << t.real() << " msecs.\n";
01159 #endif
01160 }
01161 
01162 //: computes Lucas-Kanade optical flow on the complete input views
01163 void brip_vil_float_ops::
01164 lucas_kanade_motion_on_view(vil_image_view<float> const& curr_frame,
01165                             vil_image_view<float> const& prev_frame,
01166                             const double thresh,
01167                             float& vx,
01168                             float& vy)
01169 {
01170   unsigned w = curr_frame.ni(), h = curr_frame.nj();
01171   unsigned N  = w*h;
01172   vil_image_view<float> grad_x, grad_y;
01173   grad_x.set_size(w,h);
01174   grad_y.set_size(w,h);
01175   // compute the gradient vector and the time derivative
01176   brip_vil_float_ops::gradient_3x3(curr_frame, grad_x, grad_y);
01177   vil_image_view<float> diff =
01178     brip_vil_float_ops::difference(prev_frame, curr_frame);
01179 
01180   // sum the motion terms over the view
01181   float IxIx=0.0f, IxIy=0.0f, IyIy=0.0f, IxIt=0.0f, IyIt=0.0f, dsum=0.0f;
01182   for (unsigned j = 0; j<h; j++)
01183     for (unsigned i = 0; i<w; i++)
01184     {
01185       float gx = grad_x(i, j), gy = grad_y(i, j);
01186       float dt = diff(i, j);
01187       dsum += dt*dt;
01188       IxIx += gx*gx;
01189       IxIy += gx*gy;
01190       IyIy += gy*gy;
01191       IxIt += gx*dt;
01192       IyIt += gy*dt;
01193     }
01194   // Divide by the number of pixels in the neighborhood
01195   IxIx/=(float)N;  IxIy/=(float)N; IyIy/=(float)N; IxIt/=(float)N; IyIt/=(float)N; dsum/=(float)N;
01196   float det = float(IxIx*IyIy-IxIy*IxIy);
01197   // Eliminate small motion factors
01198   float dif = vcl_sqrt(dsum);
01199   float motion_factor = vcl_fabs(det*dif);
01200   if (motion_factor<thresh)
01201   {
01202     vx = 0.0f;
01203     vy = 0.0f;
01204     return;
01205   }
01206   // solve for the motion vector
01207   vx = (IyIy*IxIt-IxIy*IyIt)/det;
01208   vy = (-IxIy*IxIt + IxIx*IyIt)/det;
01209 }
01210 
01211 //------------------------------------------------------------------------
01212 // Assume that curr match region is larger than prev_region by the required
01213 // search ranges. Step through the search and output the shift that
01214 // maximizes the correlation. zero_i and zero_j indicate the curr_image
01215 // pixel location corresponding to no velocity between frames
01216 void brip_vil_float_ops::
01217 velocity_by_correlation(vil_image_view<float> const& curr_image,
01218                         vil_image_view<float> const& prev_region,
01219                         const unsigned start_i, const unsigned end_i,
01220                         const unsigned start_j, const unsigned end_j,
01221                         const unsigned zero_i, const unsigned zero_j,
01222                         float& vx,
01223                         float& vy)
01224 {
01225   unsigned ni = prev_region.ni(), nj = prev_region.nj();
01226   float corr_max = -10.0f;
01227   float vx0 = static_cast<float>(zero_i);
01228   float vy0 = static_cast<float>(zero_j);
01229   float area = static_cast<float>(ni*nj);
01230   vx = 0; vy = 0;
01231   unsigned max_i = start_i, max_j = start_j;
01232   for (unsigned j = start_j; j<=end_j; ++j)
01233     for (unsigned i = start_i; i<=end_i; ++i)
01234     {
01235       float si1 = 0, si2 = 0, si1i1 = 0, si2i2 = 0, si1i2 = 0;
01236       //sum over the region
01237       for (unsigned r = 0; r<nj; ++r)
01238         for (unsigned c = 0; c<ni; ++c)
01239         {
01240           float I1 = prev_region(c, r);
01241           float I2 = curr_image(i+c, j+r);
01242           si1 += I1; si2 += I2;
01243           si1i1 += I1*I1;
01244           si2i2 += I2*I2;
01245           si1i2 += I1*I2;
01246         }
01247       float corr = cross_corr(area, si1, si2, si1i1, si2i2,si1i2, 1.0f);
01248       if (corr>corr_max)
01249       {
01250         corr_max = corr;
01251         max_i = i; max_j = j;
01252       }
01253 #if 0
01254       float di = i-vx0, dj = j-vy0;
01255       vcl_cout <<  di << '\t' << dj << '\t' << corr << '\n';
01256 #endif
01257     }
01258   // the velocity is given by the max indices relative to the zero location
01259   vx = static_cast<float>(max_i)- vx0;
01260   vy = static_cast<float>(max_j) - vy0;
01261   // vcl_cout << '(' << vx << ' ' << vy << "): " << corr_max << '\n';
01262 }
01263 
01264 //---------------------------------------------------------------------
01265 // Horn-Schunk method for calc. motion vectors:  Solve for the motion vectors
01266 // iteratively using a cost function with two terms (RHS of optical flow eqn
01267 // and the magnitude of spatial derivatives of the velocity field
01268 // (so that pixel-to-pixel variations are small). The second term is
01269 // weighted by alpha_coef term. The iteration goes on until the error
01270 // reaches below err_thresh
01271 //  Error conditions:
01272 //  -  -1  -  current_frame and previous_frame are equal
01273 //  -  -2  -  at least one input frame or internal process image is all zeros
01274 //  -   0  -  routine was successful
01275 
01276 int brip_vil_float_ops::
01277 Horn_SchunckMotion(vil_image_view<float> const& current_frame,
01278                    vil_image_view<float> const& previous_frame,
01279                    vil_image_view<float>& vx,
01280                    vil_image_view<float>& vy,
01281                    const float alpha_coef,
01282                    const int no_of_iterations)
01283 {
01284   // Check for equal images
01285   if (vil_image_view_deep_equality (previous_frame, current_frame ) )
01286   {
01287     vcl_cout<<"Images are same";
01288     return -1;
01289   }
01290 
01291   // Declarations
01292   unsigned w = current_frame.ni(), h = current_frame.nj();
01293 
01294   vil_image_view<float> grad_x, grad_y, diff;
01295 
01296   vil_image_view<float> temp1;
01297   vil_image_view<float> temp2;
01298 
01299   vil_image_view<float> emptyimg;
01300 
01301   // Size Init
01302 
01303   grad_x.set_size(w,h);
01304   grad_y.set_size(w,h);
01305   diff.set_size(w,h);
01306   temp1.set_size(w,h);
01307   temp2.set_size(w,h);
01308 
01309   emptyimg.set_size(w,h);
01310 
01311   temp1.fill(0.0);
01312   temp2.fill(0.0);
01313 
01314   // Initialization
01315   for (unsigned y = 0; y<h; y++)
01316     for (unsigned x = 0; x<w; x++)
01317     {
01318       vx(x,y)=0.0f;
01319       vy(x,y)=0.0f;
01320       diff (x,y)=0.0f;
01321       grad_x (x, y)= 0.0f;
01322       grad_y (x, y)= 0.0f;
01323 
01324       emptyimg (x, y) = 0.0f;
01325     }
01326 
01327   // Check for empty images
01328   if ( (vil_image_view_deep_equality (emptyimg, current_frame )) || (vil_image_view_deep_equality(emptyimg, previous_frame)))
01329   {
01330     vcl_cout<<"Image is empty";
01331     return -2;
01332   }
01333 
01334   // compute the gradient vector for current and previous
01335   brip_vil_float_ops::gradient_3x3 (current_frame , grad_x , grad_y);
01336   brip_vil_float_ops::gradient_3x3 (previous_frame , temp1 , temp2);
01337 
01338   //  Grad = 0.5* Grad(current) + 0.5 * Grad(previous)
01339   vil_math_add_image_fraction(grad_x, 0.5, temp1, 0.5);
01340   vil_math_add_image_fraction(grad_y, 0.5, temp2, 0.5);
01341   if ( (vil_image_view_deep_equality(emptyimg, grad_x)) ||
01342        (vil_image_view_deep_equality(emptyimg, grad_y)) )
01343     {
01344       vcl_cout<<"Gradient Image is empty";
01345       return -2;
01346     }
01347 
01348   temp1.fill(0.0);
01349   temp2.fill(0.0);
01350 
01351   // Average the local intensities over 3x3 region
01352   temp1 = brip_vil_float_ops::average_NxN(previous_frame, 3);
01353   temp2 = brip_vil_float_ops::average_NxN(current_frame, 3);
01354   if (vil_image_view_deep_equality(emptyimg, temp1) ||
01355       vil_image_view_deep_equality(emptyimg, temp2))
01356     {
01357       vcl_cout<<"Averaged Image is empty";
01358       return -2;
01359     }
01360 
01361   // Compute the time derivative (difference of local average intensities)
01362   // diff = dI/dt
01363   diff = brip_vil_float_ops::difference(temp1 , temp2);
01364   if (vil_image_view_deep_equality(emptyimg, diff) )
01365   {
01366     vcl_cout<<"Difference Image is empty";
01367     return -2;
01368   }
01369 
01370   temp1.fill(0.0);
01371   temp2.fill(0.0);
01372   // Iterate
01373 #ifdef DEBUG
01374   vul_timer t;
01375 #endif
01376   for (int i=0;i<no_of_iterations;i++)
01377   {
01378     // Update vx and vy
01379     //Smoothed velocities on 3x3 region
01380     temp1 = brip_vil_float_ops::average_NxN (vx,  3);
01381     temp2 = brip_vil_float_ops::average_NxN (vy,  3);
01382     for (unsigned y = 1; y+1<h; ++y)
01383       for (unsigned x = 1; x+1<w; ++x)
01384       {
01385         float tempx = temp1(x,y);
01386         float tempy = temp2(x,y);
01387 
01388         float gx = grad_x(x, y), gy = grad_y(x, y);
01389 
01390         float dt = diff(x, y);
01391         //         _____
01392         // term = (v(x,y).Grad(x,y) + dI/dt(x,y))/(alpha + |Grad(x,y)|^2)
01393         // term is the brightness constraint normalized by gradient mag.
01394         //
01395         float term =
01396           ( (gx * tempx) + (gy * tempy) + dt )/ (alpha_coef + gx*gx + gy*gy);
01397 
01398         //         ______
01399         //v(x,y) = v(x,y) - Grad(x,y)* term
01400         vx(x,y) = tempx - (gx *  term);
01401         vy(x,y) = tempy - (gy *  term);
01402       }
01403 
01404 #ifdef DEBUG
01405     vcl_cout << "Iteration No " << i << '\n';
01406 #endif
01407     brip_vil_float_ops::fill_x_border(vx, 1, 0.0f);
01408     brip_vil_float_ops::fill_y_border(vx, 1, 0.0f);
01409     brip_vil_float_ops::fill_x_border(vy, 1, 0.0f);
01410     brip_vil_float_ops::fill_y_border(vy, 1, 0.0f);
01411   }
01412 #ifdef DEBUG
01413   vcl_cout << "\nCompute Horn-Schunck iteration in " << t.real() << " msecs.\n";
01414 #endif
01415   return 0;
01416 }
01417 
01418 void brip_vil_float_ops::fill_x_border(vil_image_view<float> & image,
01419                                        unsigned w, float value)
01420 {
01421   unsigned width = image.ni(), height = image.nj();
01422   if (2*w>width)
01423   {
01424     vcl_cout << "In brip_vil_float_ops::fill_x_border(..) - 2xborder exceeds image width\n";
01425     return;
01426   }
01427   for (unsigned y = 0; y<height; y++)
01428     for (unsigned x = 0; x<w; x++)
01429       image(x, y) = value;
01430 
01431   for (unsigned y = 0; y<height; y++)
01432     for (unsigned x = width-w; x<width; x++)
01433       image(x, y) = value;
01434 }
01435 
01436 void brip_vil_float_ops::fill_y_border(vil_image_view<float> & image,
01437                                        unsigned h, float value)
01438 {
01439   unsigned width = image.ni(), height = image.nj();
01440   if (2*h>height)
01441   {
01442     vcl_cout << "In brip_vil_float_ops::fill_y_border(..) - 2xborder exceeds image height\n";
01443     return;
01444   }
01445   for (unsigned y = 0; y<h; y++)
01446     for (unsigned x = 0; x<width; x++)
01447       image(x, y) = value;
01448 
01449   for (unsigned y = height-h; y<height; y++)
01450     for (unsigned x = 0; x<width; x++)
01451       image(x, y) = value;
01452 }
01453 
01454 vil_image_view<vxl_byte>
01455 brip_vil_float_ops::convert_to_byte(vil_image_view<float> const& image)
01456 {
01457   // determine the max min values
01458   float min_val = vnl_numeric_traits<float>::maxval;
01459   float max_val = -min_val;
01460   unsigned w = image.ni(), h = image.nj();
01461   vil_image_view<vxl_byte> output;
01462   output.set_size(w,h);
01463   for (unsigned y = 0; y<h; y++)
01464     for (unsigned x = 0; x<w; x++)
01465     {
01466       min_val = vnl_math::min(min_val, image(x,y));
01467       max_val = vnl_math::max(max_val, image(x,y));
01468     }
01469   float range = max_val-min_val;
01470   if (range == 0.f)
01471     range = 1.f;
01472   else
01473     range = 255.f/range;
01474   for (unsigned y = 0; y<h; y++)
01475     for (unsigned x = 0; x<w; x++)
01476     {
01477       float v = (image(x,y)-min_val)*range;
01478       if(v>255.0f) v = 255.0f;//could have small round off above 255
01479       output(x,y) = (unsigned char)v;
01480     }
01481   return output;
01482 }
01483 
01484 //------------------------------------------------------------
01485 //: Convert the range between min_val and max_val to 255
01486 vil_image_view<vxl_byte>
01487 brip_vil_float_ops::convert_to_byte(vil_image_view<float> const& image,
01488                                     float min_val, float max_val)
01489 {
01490   unsigned w = image.ni(), h = image.nj();
01491   vil_image_view<vxl_byte> output;
01492   output.set_size(w,h);
01493   float range = max_val-min_val;
01494   if (range == 0.f)
01495     range = 1.f;
01496   else
01497     range = 256.f/range;
01498   for (unsigned y = 0; y<h; y++)
01499     for (unsigned x = 0; x<w; x++)
01500     {
01501       float v = (image(x,y)-min_val)*range;
01502       if (v>=256)
01503         v=255;
01504       else if (v<0)
01505         v=0;
01506       output(x,y) = (unsigned char)v;
01507     }
01508 
01509   return output;
01510 }
01511 
01512 vil_image_view<vxl_byte> brip_vil_float_ops::
01513 convert_to_byte(vil_image_view<vxl_uint_16> const& image,
01514                 vxl_uint_16 min_val, vxl_uint_16 max_val)
01515 {
01516   unsigned ni = image.ni(), nj = image.nj();
01517   vil_image_view<vxl_byte> output;
01518   output.set_size(ni, nj);
01519   float range = static_cast<float>(max_val-min_val);
01520   if (!range)
01521     range = 1.f;
01522   else
01523     range = 256.f/range;
01524   for (unsigned r = 0; r<nj; r++)
01525     for (unsigned c = 0; c<ni; c++)
01526     {
01527       float v = float(image(c, r)-min_val)*range;
01528       if (v>=256)
01529         v=255;
01530       else if (v<0)
01531         v=0;
01532       output(c, r) = static_cast<unsigned char>(v);
01533     }
01534   return output;
01535 }
01536 
01537 // Note this is a more standard interface than convert_to_grey
01538 vil_image_view<vxl_byte>
01539 brip_vil_float_ops::convert_to_byte(vil_image_resource_sptr const& image)
01540 {
01541   return brip_vil_float_ops::convert_to_grey(*image);
01542 }
01543 
01544 vil_image_view<vxl_uint_16>
01545 brip_vil_float_ops::convert_to_short(vil_image_view<float> const& image,
01546                                      float min_val, float max_val)
01547 {
01548   unsigned w = image.ni(), h = image.nj();
01549   float max_short = 65355.f;
01550   vil_image_view<vxl_uint_16> output;
01551   output.set_size(w,h);
01552   float range = max_val-min_val;
01553   if (!range)
01554     range = 1;
01555   else
01556     range = max_short/range;
01557   for (unsigned y = 0; y<h; y++)
01558     for (unsigned x = 0; x<w; x++)
01559     {
01560       float v = (image(x,y)-min_val)*range;
01561       if (v>max_short)
01562         v=max_short;
01563       if (v<0)
01564         v=0;
01565       output(x,y) = (unsigned short)v;
01566     }
01567   return output;
01568 }
01569 
01570 //: converts a float image to an unsigned short (16-bit) image.
01571 // range determined automatically
01572 vil_image_view<vxl_uint_16>
01573 brip_vil_float_ops::convert_to_short(vil_image_view<float> const& image)
01574 {
01575   float minv = vnl_numeric_traits<float>::maxval;
01576   float maxv = -minv;
01577   unsigned ni = image.ni(), nj = image.nj();
01578   for (unsigned j = 0; j<nj; ++j)
01579     for (unsigned i = 0; i<ni; ++i)
01580     {
01581       float v = image(i,j);
01582       if (v<minv)
01583         minv = v;
01584       if (v>maxv)
01585         maxv = v;
01586     }
01587   return brip_vil_float_ops::convert_to_short(image, minv, maxv);
01588 }
01589 
01590 vil_image_view<vxl_uint_16>
01591 brip_vil_float_ops::convert_to_short(vil_image_resource_sptr const& image)
01592 {
01593   // Check if the image is a float
01594   if (image->nplanes()==1 &&image->pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
01595   {
01596     vil_image_view<float> temp = image->get_view();
01597     float vmin=0, vmax= 65355;
01598     vil_math_value_range<float>(temp, vmin, vmax);
01599     return brip_vil_float_ops::convert_to_short(temp, vmin, vmax);
01600   }
01601 
01602   // Here we assume that the image is an unsigned char
01603   if (image->nplanes()==1&&image->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
01604   {
01605     vil_image_view<unsigned char > temp = image->get_view();
01606     vil_image_view<unsigned short> short_image;
01607     unsigned width = temp.ni(), height = temp.nj();
01608     short_image.set_size(width, height);
01609     for (unsigned y = 0; y<height; y++)
01610       for (unsigned x = 0; x<width; x++)
01611         short_image(x,y) = static_cast<unsigned short>(temp(x,y));
01612     return temp;
01613   }
01614 
01615   // Here the image is an unsigned short (16-bit) image so just return it
01616   if (image->nplanes()==1&&image->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
01617   {
01618     vil_image_view<unsigned short > temp = image->get_view();
01619     return temp;
01620   }
01621 
01622   // the image is color so we should convert it to greyscale
01623   // Here we assume the color elements are unsigned char.
01624   if (image->nplanes()==3&&image->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
01625   {
01626     vil_image_view<vxl_byte> color_image = image->get_view();
01627     unsigned width = color_image.ni(), height = color_image.nj();
01628     // the output image
01629     vil_image_view<unsigned short> short_image;
01630     short_image.set_size(width, height);
01631     for (unsigned y = 0; y<height; y++)
01632       for (unsigned x = 0; x<width; x++)
01633       {
01634         double v = color_image(x,y,0)+color_image(x,y,1)+color_image(x,y,2);
01635         v/=3.0;
01636         short_image(x,y) = static_cast<unsigned short>(v);
01637       }
01638     return short_image;
01639   }
01640 
01641   // the image is multispectral so we should convert it to greyscale
01642   // Here we assume the color elements are unsigned short (16-bit).
01643   if (image->nplanes()==4&&image->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
01644   {
01645     vil_image_view<unsigned short > mband_image = image->get_view();
01646     unsigned width = mband_image.ni(), height = mband_image.nj();
01647     // the output image
01648     vil_image_view<unsigned short> short_image;
01649     short_image.set_size(width, height);
01650     for (unsigned y = 0; y<height; y++)
01651       for (unsigned x = 0; x<width; x++)
01652       {
01653         unsigned short v = 0;
01654         for (unsigned p = 0; p<4; ++p)
01655           v += mband_image(x, y, p);
01656         v/=4;
01657         short_image(x,y) = v;
01658       }
01659     return short_image;
01660   }
01661 
01662   // If we get here then the input is not a type we handle so return a null view
01663   return vil_image_view<vxl_uint_16>();
01664 }
01665 
01666 vil_image_view<float>
01667 brip_vil_float_ops::convert_to_float(vil_image_view<vxl_byte> const& image)
01668 {
01669   vil_image_view<float> output;
01670   unsigned ni = image.ni(), nj = image.nj(), np = image.nplanes();
01671   output.set_size(ni,nj, np);
01672   for (unsigned j = 0; j<nj; j++)
01673     for (unsigned i = 0; i<ni; i++)
01674       for (unsigned p = 0; p<np; ++p)
01675         output(i,j,p) = (float)image(i,j,p);
01676   return output;
01677 }
01678 
01679 vil_image_view<float>
01680 brip_vil_float_ops::convert_to_float(vil_image_view<vxl_uint_16> const& image)
01681 {
01682   vil_image_view<float> output;
01683   unsigned w = image.ni(), h = image.nj();
01684   output.set_size(w,h);
01685   for (unsigned y = 0; y<h; y++)
01686     for (unsigned x = 0; x<w; x++)
01687       output(x,y) = (float)image(x,y);
01688   return output;
01689 }
01690 
01691 vil_image_view<bool>
01692 brip_vil_float_ops::convert_to_bool(vil_image_view<vxl_byte> const& image)
01693 {
01694   vil_image_view<bool> output;
01695   unsigned w = image.ni(), h = image.nj();
01696   output.set_size(w,h);
01697   for (unsigned y = 0; y<h; y++)
01698     for (unsigned x = 0; x<w; x++)
01699       if (image(x,y) > 128)
01700         output(x,y)=true;
01701       else
01702         output(x,y)=false;
01703   return output;
01704 }
01705 
01706 vil_image_view<float>
01707 brip_vil_float_ops::convert_to_float(vil_image_view<vil_rgb<vxl_byte> > const& image)
01708 {
01709   vil_image_view<float> output;
01710   unsigned w = image.ni(), h = image.nj();
01711   output.set_size(w,h);
01712   for (unsigned y = 0; y<h; y++)
01713     for (unsigned x = 0; x<w; x++)
01714     {
01715       vil_rgb<vxl_byte> rgb = image(x,y);
01716       output(x,y) = (float)rgb.grey();
01717     }
01718   return output;
01719 }
01720 
01721 void brip_vil_float_ops::rgb_to_ihs(vil_rgb<vxl_byte> const& rgb,
01722                                     float& i, float& h, float& s)
01723 {
01724   float r=rgb.R();
01725   float g=rgb.G();
01726   float b=rgb.B();
01727 
01728   float maxval = vnl_math::max(r,vnl_math::max(g,b));
01729   float minval = vnl_math::min(r,vnl_math::min(g,b));
01730 
01731   float delta = maxval - minval;
01732   i = maxval;
01733   if (maxval == 0)
01734     s = 0;
01735   else
01736     s = delta / maxval;
01737 
01738   if (s==0)
01739     h = 0;                   //!< (Hue is undefined)
01740 
01741   if (r== maxval)
01742     h = (g - b) / delta ;    //!< (between yellow and magenta)
01743   if (g == maxval)
01744     h = 2 + (b - r)/delta ;  //!< (between cyan and yellow)
01745   if (b == maxval)
01746     h = 4 + (r - g) / delta; //!< (between magenta and cyan)
01747   h *= 60.f;                 //!< (convert Hue to degrees)
01748   if (h < 0.f)
01749     h += 360.f;              //!< (Hue must be positive)
01750   else if (h >= 360.f)
01751     h -= 360.f;              //!< (Hue must be less than 360)
01752 
01753   h *= 256.0f / 360.0f; // range 0 .. 255
01754   s *= 256.0f;          // range 0 .. 255
01755 }
01756 
01757 void brip_vil_float_ops::rgb_to_ihs_tsai(vil_rgb<vxl_byte> const& rgb, float& i, float& h, float& s)
01758 {
01759   float r=rgb.R();
01760   float g=rgb.G();
01761   float b=rgb.B();
01762   float sq6 = vcl_sqrt(6.0f);
01763   float V1, V2;
01764   i = (r+g+b)/3.0f;
01765   V1 = (2.0f*b-r-g)/sq6;
01766   V2 = (r-2.0f*g)/sq6;
01767   s = vcl_sqrt(V1*V1 + V2*V2);
01768   float two_pi = static_cast<float>(vnl_math::twopi);
01769   float a = vcl_atan2(V2, V1)/two_pi;
01770   if (a<0.0f) a += 1.0f; // range [0..1)
01771   h = a*256.0f; // range [0..256)
01772 }
01773 
01774 void brip_vil_float_ops::ihs_to_rgb(vil_rgb<vxl_byte> & rgb,
01775                                     const float i, const float h, const float s)
01776 {
01777   // Reference: page 593 of Foley & van Dam
01778   float R = 0.0f;
01779   float G = 0.0f;
01780   float B = 0.0f;
01781 
01782   if (s == 0) {
01783     R=i;
01784     G=i;
01785     B=i;
01786   }
01787   else if (s > 0.0)
01788   {
01789     float ss = s, hh = h;
01790     ss *= 1.f / 256.f;
01791     hh *= 6.f / 256.f;
01792 
01793     float J = vcl_floor(hh);
01794     float F = hh - J;
01795     float P =( i * (1 - ss));
01796     float Q = (i * (1 - (ss * F)));
01797     float T = (i * (1 - (ss * (1 - F))));
01798 
01799     if (J == 0) { R=i; G=T; B=P; }
01800     if (J == 1) { R=Q; G=i; B=P; }
01801     if (J == 2) { R=P; G=i; B=T; }
01802     if (J == 3) { R=P; G=Q; B=i; }
01803     if (J == 4) { R=T; G=P; B=i; }
01804     if (J == 5) { R=i; G=P; B=Q; }
01805   }
01806   rgb.r = (vxl_byte)R;
01807   rgb.g = (vxl_byte)G;
01808   rgb.b = (vxl_byte)B;
01809 }
01810 
01811 void brip_vil_float_ops::
01812 convert_to_IHS(vil_image_view<vil_rgb<vxl_byte> > const& image,
01813                vil_image_view<float>& I,
01814                vil_image_view<float>& H,
01815                vil_image_view<float>& S)
01816 {
01817   unsigned w = image.ni(), h = image.nj();
01818   I.set_size(w,h);
01819   H.set_size(w,h);
01820   S.set_size(w,h);
01821   for (unsigned r = 0; r < h; r++)
01822     for (unsigned c = 0; c < w; c++)
01823     {
01824       float in, hue, sat;
01825       rgb_to_ihs(image(c,r), in, hue, sat);
01826       I(c,r) = in;
01827       H(c,r) = hue;
01828       S(c,r) = sat;
01829     }
01830 }
01831 
01832 void brip_vil_float_ops::
01833 convert_to_IHS(vil_image_view<vxl_byte> const& image,
01834                vil_image_view<float>& I,
01835                vil_image_view<float>& H,
01836                vil_image_view<float>& S)
01837 {
01838   unsigned w = image.ni(), h = image.nj();
01839   I.set_size(w,h);
01840   H.set_size(w,h);
01841   S.set_size(w,h);
01842   for (unsigned r = 0; r < h; r++)
01843     for (unsigned c = 0; c < w; c++)
01844     {
01845       float in, hue, sat;
01846       vil_rgb<vxl_byte> imint(image(c,r,0),image(c,r,1),image(c,r,2));
01847       rgb_to_ihs(imint, in, hue, sat);
01848       I(c,r) = in;
01849       H(c,r) = hue;
01850       S(c,r) = sat;
01851     }
01852 }
01853 
01854 void brip_vil_float_ops::
01855 convert_to_IHS_tsai(vil_image_view<vxl_byte> const& image,
01856                     vil_image_view<float>& I,
01857                     vil_image_view<float>& H,
01858                     vil_image_view<float>& S,
01859                     vil_image_view<float>& ratio)
01860 {
01861   unsigned w = image.ni(), h = image.nj();
01862   I.set_size(w,h);
01863   H.set_size(w,h);
01864   S.set_size(w,h);
01865   ratio.set_size(w,h);
01866   for (unsigned r = 0; r < h; r++)
01867     for (unsigned c = 0; c < w; c++)
01868     {
01869       float in, hue, sat;
01870       vil_rgb<vxl_byte> imint(image(c,r,0),image(c,r,1),image(c,r,2));
01871       rgb_to_ihs_tsai(imint, in, hue, sat);
01872       I(c,r) = in;
01873       H(c,r) = hue;
01874       S(c,r) = sat;
01875       ratio(c,r) = (hue+1.0f)/(in+1.0f);
01876     }
01877 }
01878 
01879 #if 0 // commented out
01880 void brip_vil_float_ops::
01881 display_IHS_as_RGB(vil_image_view<float> const& I,
01882                    vil_image_view<float> const& H,
01883                    vil_image_view<float> const& S,
01884                    vil_image_view<vil_rgb<vxl_byte> >& image)
01885 {
01886   unsigned w = I.ni(), h = I.nj();
01887   image.set_size(w,h);
01888   float s = 256.0f/360.0f;
01889   for (unsigned r = 0; r < h; r++)
01890     for (unsigned c = 0; c < w; c++)
01891     {
01892       float in = I(c,r);
01893       float hue = s * H(c,r);
01894       float sat = S(c,r);
01895       if (in<0) in = 0;
01896       if (sat<0) sat = 0;
01897       if (hue<0) hue = 0;
01898       if (in>=256) in = 255;
01899       if (hue>=256) hue = 255;
01900       if (sat>=256) sat = 255;
01901       image(c,r).r = (vxl_byte)in;
01902       image(c,r).g = (vxl_byte)hue;
01903       image(c,r).b = (vxl_byte)sat;
01904     }
01905 }
01906 #endif // 0
01907 
01908 // map so that intensity is proportional to saturation and hue is color
01909 void brip_vil_float_ops::
01910 display_IHS_as_RGB(vil_image_view<float> const& I,
01911                    vil_image_view<float> const& H,
01912                    vil_image_view<float> const& S,
01913                    vil_image_view<vil_rgb<vxl_byte> >& image)
01914 {
01915   unsigned w = I.ni(), h = I.nj();
01916   image.set_size(w,h);
01917 
01918   const float deg_to_rad = float(vnl_math::pi_over_180);
01919   for (unsigned r = 0; r < h; r++)
01920     for (unsigned c = 0; c < w; c++)
01921     {
01922       float hue = H(c,r);
01923       float sat = 2.f*S(c,r);
01924       if (sat<0)
01925         sat = 0.f;
01926       else if (sat>=256)
01927         sat = 255.999f;
01928       float ang = deg_to_rad*hue;
01929       float cs = vcl_cos(ang), si = vcl_fabs(vcl_sin(ang));
01930       float red=0.0f, blue=0.0f;
01931       float green = si*sat;
01932       if (cs>=0)
01933         red = cs*sat;
01934       else
01935         blue = sat*(-cs);
01936       image(c,r).r = (vxl_byte)red;
01937       image(c,r).g = (vxl_byte)green;
01938       image(c,r).b = (vxl_byte)blue;
01939     }
01940 }
01941 
01942 vil_image_view<vil_rgb<vxl_byte> > brip_vil_float_ops::
01943 combine_color_planes(vil_image_view<vxl_byte> const& R,
01944                      vil_image_view<vxl_byte> const& G,
01945                      vil_image_view<vxl_byte> const& B)
01946 {
01947   unsigned w = R.ni(), h = R.nj();
01948   vil_image_view<vil_rgb<vxl_byte> > image(w,h);
01949   for (unsigned r = 0; r < h; r++)
01950     for (unsigned c = 0; c < w; c++)
01951     {
01952       image(c,r).r = R(c,r);
01953       image(c,r).g = G(c,r);
01954       image(c,r).b = B(c,r);
01955     }
01956   return image;
01957 }
01958 
01959 vil_image_view<vil_rgb<vxl_byte> >
01960 brip_vil_float_ops::combine_color_planes(vil_image_resource_sptr const& R,
01961                                          vil_image_resource_sptr const& G,
01962                                          vil_image_resource_sptr const& B)
01963 {
01964   vil_image_view<vil_rgb<vxl_byte> > view(0,0);
01965   if (!R||!G||!B)
01966     return view; // return an empty view
01967   // determine the union of all the resources
01968   vbl_bounding_box<unsigned int, 2> b;
01969   unsigned int r_ni = R->ni(), r_nj = R->nj();
01970   unsigned int g_ni = G->ni(), g_nj = G->nj();
01971   unsigned int b_ni = B->ni(), b_nj = B->nj();
01972   b.update(r_ni, r_nj);
01973   b.update(g_ni, g_nj);
01974   b.update(b_ni, b_nj);
01975   unsigned int n_i = b.xmax(), n_j = b.ymax();
01976   view.set_size(n_i, n_j);
01977   vil_rgb<vxl_byte> zero(0, 0, 0);
01978   vil_image_view<float>    fR = brip_vil_float_ops::convert_to_float(R);
01979   vil_image_view<vxl_byte> cR = brip_vil_float_ops::convert_to_byte(fR);
01980   vil_image_view<float>    fG = brip_vil_float_ops::convert_to_float(G);
01981   vil_image_view<vxl_byte> cG = brip_vil_float_ops::convert_to_byte(fG);
01982   vil_image_view<float>    fB = brip_vil_float_ops::convert_to_float(B);
01983   vil_image_view<vxl_byte> cB = brip_vil_float_ops::convert_to_byte(fB);
01984   for (vxl_uint_16 j = 0; j<n_j; ++j)
01985     for (vxl_uint_16 i = 0; i<n_i; ++i)
01986     {
01987       vil_rgb<vxl_byte> v = zero;
01988       if (i<r_ni&&j<r_nj)
01989         v.r = cR(i,j);
01990       if (i<g_ni&&j<g_nj)
01991         v.g= cG(i,j);
01992       if (i<b_ni&&j<b_nj)
01993         v.b= cB(i,j);
01994       view(i,j)=v;
01995     }
01996   return view;
01997 }
01998 
01999 vil_image_view<float>
02000 brip_vil_float_ops::convert_to_float(vil_image_resource const& image)
02001 {
02002   vil_image_view<float> fimg;
02003   vil_pixel_format fmt = image.pixel_format();
02004   if (vil_pixel_format_num_components(image.pixel_format())==1)
02005   {
02006     if (fmt==VIL_PIXEL_FORMAT_UINT_16)
02007     {
02008       vil_image_view<unsigned short> temp=image.get_view();
02009       fimg = brip_vil_float_ops::convert_to_float(temp);
02010     }
02011     else if (fmt==VIL_PIXEL_FORMAT_BYTE)
02012     {
02013       vil_image_view<unsigned char> temp=image.get_view();
02014       fimg = brip_vil_float_ops::convert_to_float(temp);
02015     }
02016     else if (fmt==VIL_PIXEL_FORMAT_FLOAT)
02017       return image.get_view();
02018   }
02019   else if (vil_pixel_format_num_components(fmt)==3)
02020   {
02021     vil_image_view<vil_rgb<vxl_byte> > temp= image.get_view();
02022     fimg = brip_vil_float_ops::convert_to_float(temp);
02023   }
02024   else
02025   {
02026     vcl_cout << "In brip_vil_float_ops::convert_to_float - input not color or grey\n";
02027     return vil_image_view<float>();
02028   }
02029   return fimg;
02030 }
02031 
02032 //-----------------------------------------------------------------
02033 //: Convert any image to an unsigned_char image
02034 vil_image_view<vxl_byte>
02035 brip_vil_float_ops::convert_to_grey(vil_image_resource const& image)
02036 {
02037   // Check if the image is a float
02038   if (image.nplanes()==1 &&image.pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
02039   {
02040     vil_image_view<float> temp = image.get_view();
02041     float vmin=0.f, vmax=256.f;
02042     vil_math_value_range<float>(temp, vmin, vmax);
02043     return brip_vil_float_ops::convert_to_byte(temp, vmin, vmax);
02044   }
02045   if (image.nplanes()==1 &&image.pixel_format()==VIL_PIXEL_FORMAT_BOOL)
02046   {
02047     vil_image_view<bool> temp = image.get_view();
02048     unsigned nj = temp.nj(), ni = temp.ni();
02049     vil_image_view<unsigned char> out(ni, nj);
02050     out.fill(0);
02051     for (unsigned j = 0; j<nj; ++j)
02052       for (unsigned i = 0; i<ni; ++i)
02053         if (temp(i,j)) out(i,j) = 255;
02054     return out;
02055   }
02056 
02057   // Here we assume that the image is an unsigned char
02058   // In this case we should just return it.
02059   if (image.nplanes()==1&&image.pixel_format()==VIL_PIXEL_FORMAT_BYTE)
02060   {
02061     vil_image_view<unsigned char > temp = image.get_view();
02062     return temp;
02063   }
02064 
02065   if (image.nplanes()==1&&image.pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
02066   {
02067     vil_image_view<unsigned short > temp = image.get_view();
02068     unsigned short vmin=0, vmax=255;
02069     vil_math_value_range<unsigned short>(temp, vmin, vmax);
02070     return brip_vil_float_ops::convert_to_byte(temp, vmin, vmax);
02071   }
02072   // the image is color so we should convert it to greyscale
02073   // Here we assume the color elements are unsigned char.
02074   if (image.nplanes()==3&&image.pixel_format()==VIL_PIXEL_FORMAT_BYTE)
02075   {
02076     vil_image_view<vil_rgb<vxl_byte> > color_image = image.get_view();
02077     unsigned width = color_image.ni(), height = color_image.nj();
02078     // the output image
02079     vil_image_view<unsigned char> grey_image;
02080     grey_image.set_size(width, height);
02081     for (unsigned y = 0; y<height; y++)
02082       for (unsigned x = 0; x<width; x++)
02083         grey_image(x,y) = color_image(x,y).grey();
02084     return grey_image;
02085   }
02086   if (image.nplanes()==3&&image.pixel_format()==VIL_PIXEL_FORMAT_FLOAT)
02087   {
02088     vil_image_view<float> color_image = image.get_view();
02089     unsigned width = color_image.ni(), height = color_image.nj();
02090     // the output image
02091     vil_image_view<float> grey_image_f;
02092     grey_image_f.set_size(width, height);
02093     for (unsigned y = 0; y<height; y++)
02094       for (unsigned x = 0; x<width; x++) {
02095         float v = 0;
02096         for (unsigned p = 0; p<3; ++p)
02097           v += color_image(x,y,p);
02098         grey_image_f(x,y) = v/3.0f;
02099       }
02100     float vmin=0.f, vmax=256.f;
02101     vil_math_value_range<float>(grey_image_f, vmin, vmax);
02102     return brip_vil_float_ops::convert_to_byte(grey_image_f, vmin, vmax);
02103   }
02104   // If we get here then the input is not a type we handle so return a null view
02105   return vil_image_view<vxl_byte>();
02106 }
02107 
02108 //--------------------------------------------------------------
02109 //: Read a convolution kernel from file
02110 // Assumes a square kernel with odd dimensions, i.e., w,h = 2n+1
02111 // format:
02112 // \verbatim
02113 //     n
02114 //     scale
02115 //     k00  k01  ... k02n
02116 //           ...
02117 //     k2n0 k2n1 ... k2n2n
02118 // \endverbatim
02119 //
02120 vbl_array_2d<float> brip_vil_float_ops::load_kernel(vcl_string const& file)
02121 {
02122   vcl_ifstream instr(file.c_str(), vcl_ios::in);
02123   if (!instr)
02124   {
02125     vcl_cout << "In brip_vil_float_ops::load_kernel - failed to load kernel\n";
02126     return vbl_array_2d<float>(0,0);
02127   }
02128   unsigned n;
02129   float scale;
02130   float v =0.0f;
02131   instr >> n;
02132   instr >> scale;
02133   unsigned N = 2*n+1;
02134   vbl_array_2d<float> output(N, N);
02135   for (unsigned y = 0; y<N; y++)
02136     for (unsigned x = 0; x<N; x++)
02137     {
02138       instr >> v;
02139       output.put(x, y, v/scale);
02140     }
02141 #ifdef DEBUG
02142   vcl_cout << "The Kernel\n";
02143   for (unsigned y = 0; y<N; y++)
02144   {
02145     for (unsigned x = 0; x<N; x++)
02146       vcl_cout << ' ' <<  output[x][y];
02147     vcl_cout << '\n';
02148   }
02149 #endif
02150   return output;
02151 }
02152 
02153 static void insert_image(vil_image_view<float> const& image, int col,
02154                          vnl_matrix<float> & I)
02155 {
02156   unsigned width = image.ni(), height = image.nj(), row=0;
02157   for (unsigned y =0; y<height; y++)
02158     for (unsigned x = 0; x<width; x++, row++)
02159       I.put(row, col, image(x,y));
02160 }
02161 
02162 void brip_vil_float_ops::
02163 basis_images(vcl_vector<vil_image_view<float> > const& input_images,
02164              vcl_vector<vil_image_view<float> >& basis)
02165 {
02166   basis.clear();
02167   unsigned n_images = input_images.size();
02168   if (!n_images)
02169   {
02170     vcl_cout << "In brip_vil_float_ops::basis_images(.) - no input images\n";
02171     return;
02172   }
02173   unsigned width = input_images[0].ni(), height = input_images[0].nj();
02174   unsigned npix = width*height;
02175 
02176   // Insert the images into matrix I
02177   vnl_matrix<float> I(npix, n_images, 0.f);
02178   for (unsigned i = 0; i<n_images; i++)
02179     insert_image(input_images[i], i, I);
02180 
02181   // Compute the SVD of matrix I
02182 #ifdef DEBUG
02183   vcl_cout << "Computing Singular values of a " <<  npix << " by "
02184            << n_images << " matrix\n";
02185   vul_timer t;
02186 #endif
02187   vnl_svd<float> svd(I);
02188 #ifdef DEBUG
02189   vcl_cout << "SVD Took " << t.real() << " msecs\n"
02190            << "Eigenvalues:\n";
02191   for (unsigned i = 0; i<n_images; i++)
02192     vcl_cout << svd.W(i) << '\n';
02193 #endif
02194   // Extract the Basis images
02195   unsigned rank = svd.rank();
02196   if (!rank)
02197   {
02198     vcl_cout << "In brip_vil_float_ops::basis_images(.) - I has zero rank\n";
02199     return;
02200   }
02201   vnl_matrix<float> U = svd.U();
02202   // Output the basis images
02203   unsigned rows = U.rows();
02204   for (unsigned k = 0; k<rank; k++)
02205   {
02206     vil_image_view<float> out(width, height);
02207     unsigned x =0, y = 0;
02208     for (unsigned r = 0; r<rows; r++)
02209     {
02210       out(x, y) = U(r,k);
02211       x++;
02212       if (x>=width)
02213       {
02214         y++;
02215         x=0;
02216       }
02217       if (y>=width)
02218       {
02219         vcl_cout<<"In brip_vil_float_ops::basis_images(.) - shouldn't happen\n";
02220         return;
02221       }
02222     }
02223     basis.push_back(out);
02224   }
02225 }
02226 
02227 //: 1d fourier transform
02228 //-------------------------------------------------------------------------
02229 // This computes an in-place complex-to-complex FFT
02230 // x and y are the real and imaginary arrays of 2^m points.
02231 // dir =  1 gives forward transform
02232 // dir = -1 gives reverse transform
02233 //
02234 //   Formula: forward
02235 // \verbatim
02236 //                N-1
02237 //                ---
02238 //            1   \          - j k 2 pi n / N
02239 //    X(n) = ---   >   x(k) e                    = forward transform
02240 //            N   /                                n=0..N-1
02241 //                ---
02242 //                k=0
02243 // \endverbatim
02244 //
02245 //    Formula: reverse
02246 // \verbatim
02247 //                N-1
02248 //                ---
02249 //                \          j k 2 pi n / N
02250 //    X(n) =       >   x(k) e                    = forward transform
02251 //                /                                n=0..N-1
02252 //                ---
02253 //                k=0
02254 // \endverbatim
02255 //
02256 bool brip_vil_float_ops::fft_1d(int dir, int m, double* x, double* y)
02257 {
02258   long nn,i,i1,j,k,i2,l,l1,l2;
02259   double c1,c2,tx,ty,t1,t2,u1,u2,z;
02260 
02261   // Calculate the number of points
02262   nn = 1;
02263   for (i=0;i<m;i++)
02264     nn *= 2;
02265 
02266   // Do the bit reversal
02267   i2 = nn >> 1;
02268   j = 0;
02269   for (i=0;i<nn-1;i++) {
02270     if (i < j) {
02271       tx = x[i];
02272       ty = y[i];
02273       x[i] = x[j];
02274       y[i] = y[j];
02275       x[j] = tx;
02276       y[j] = ty;
02277     }
02278     k = i2;
02279     while (k <= j) {
02280       j -= k;
02281       k >>= 1;
02282     }
02283     j += k;
02284   }
02285 
02286   // Compute the FFT
02287   c1 = -1.0;
02288   c2 = 0.0;
02289   l2 = 1;
02290   for (l=0;l<m;l++) {
02291     l1 = l2;
02292     l2 <<= 1;
02293     u1 = 1.0;
02294     u2 = 0.0;
02295     for (j=0;j<l1;j++) {
02296       for (i=j;i<nn;i+=l2) {
02297         i1 = i + l1;
02298         t1 = u1 * x[i1] - u2 * y[i1];
02299         t2 = u1 * y[i1] + u2 * x[i1];
02300         x[i1] = x[i] - t1;
02301         y[i1] = y[i] - t2;
02302         x[i] += t1;
02303         y[i] += t2;
02304       }
02305       z =  u1 * c1 - u2 * c2;
02306       u2 = u1 * c2 + u2 * c1;
02307       u1 = z;
02308     }
02309     c2 = vcl_sqrt((1.0 - c1) / 2.0);
02310     if (dir == 1)
02311       c2 = -c2;
02312     c1 = vcl_sqrt((1.0 + c1) / 2.0);
02313   }
02314 
02315   // Scaling for forward transform
02316   if (dir == 1) {
02317     for (i=0;i<nn;i++) {
02318       x[i] /= (double)nn;
02319       y[i] /= (double)nn;
02320     }
02321   }
02322 
02323   return true;
02324 }
02325 
02326 //-------------------------------------------------------------------------
02327 //: Perform a 2D FFT inplace given a complex 2D array
02328 //  The direction dir, 1 for forward, -1 for reverse
02329 //  The size of the array (nx,ny)
02330 //  Return false if there are memory problems or
02331 //  the dimensions are not powers of 2
02332 //
02333 bool brip_vil_float_ops::fft_2d(vnl_matrix<vcl_complex<double> >& c,int nx,int ny,int dir)
02334 {
02335   int i,j;
02336   int mx, my;
02337   double *real,*imag;
02338   vnl_fft_prime_factors<double> pfx (nx);
02339   vnl_fft_prime_factors<double> pfy (ny);
02340   mx = (int)pfx.pqr()[0];
02341   my = (int)pfy.pqr()[0];
02342   // Transform the rows
02343   real = new double[nx];
02344   imag = new double[nx];
02345   if (real == 0 || imag == 0)
02346     return false;
02347   for (j=0;j<ny;j++) {
02348     for (i=0;i<nx;i++) {
02349       real[i] = c[j][i].real();
02350       imag[i] = c[j][i].imag();
02351     }
02352     brip_vil_float_ops::fft_1d(dir,mx,real,imag);
02353     for (i=0;i<nx;i++) {
02354       vcl_complex<double> v(real[i], imag[i]);
02355       c[j][i] = v;
02356     }
02357   }
02358   delete [] real;
02359   delete [] imag;
02360   // Transform the columns
02361   real = new double[ny];
02362   imag = new double[ny];
02363   if (real == 0 || imag == 0)
02364     return false;
02365   for (i=0;i<nx;i++) {
02366     for (j=0;j<ny;j++) {
02367       real[j] = c[j][i].real();
02368       imag[j] = c[j][i].imag();
02369     }
02370     fft_1d(dir,my,real,imag);
02371     for (j=0;j<ny;j++) {
02372       vcl_complex<double> v(real[j], imag[j]);
02373       c[j][i] =  v;
02374     }
02375   }
02376   delete [] real;
02377   delete [] imag;
02378   return true;
02379 }
02380 
02381 //: reorder the transform values to sequential frequencies as in conventional Fourier transforms.
02382 //  The transformation is its self-inverse.
02383 void brip_vil_float_ops::
02384 ftt_fourier_2d_reorder(vnl_matrix<vcl_complex<double> > const& F1,
02385                        vnl_matrix<vcl_complex<double> > & F2)
02386 {
02387   int rows = static_cast<int>(F1.rows()), cols = static_cast<int>(F1.cols());
02388   int half_rows = rows/2, half_cols = cols/2;
02389   int ri, ci;
02390   for (int r = 0; r<rows; r++)
02391   {
02392     if (r<half_rows)
02393       ri = half_rows+r;
02394     else
02395       ri = r-half_rows;
02396     for (int c = 0; c<cols; c++)
02397     {
02398       if (c<half_cols)
02399         ci = half_cols+c;
02400       else
02401         ci = c-half_cols;
02402       F2[ri][ci]=F1[r][c];
02403     }
02404   }
02405 }
02406 
02407 //: Compute the fourier transform.
02408 //  If the image dimensions are not a power of 2 then the operation fails.
02409 bool brip_vil_float_ops::
02410 fourier_transform(vil_image_view<float> const& input,
02411                   vil_image_view<float>& mag,
02412                   vil_image_view<float>& phase)
02413 {
02414   unsigned w = input.ni(), h = input.nj();
02415 
02416   vnl_fft_prime_factors<float> pfx (w);
02417   vnl_fft_prime_factors<float> pfy (h);
02418   if (!pfx.pqr()[0]||!pfy.pqr()[0])
02419     return false;
02420   // fill the fft matrix
02421   vnl_matrix<vcl_complex<double> > fft_matrix(h, w), fourier_matrix(h,w);
02422   for (unsigned y = 0; y<h; y++)
02423     for (unsigned x =0; x<w; x++)
02424     {
02425       vcl_complex<double> cv(input(x,y), 0.0);
02426       fft_matrix.put(y, x, cv);
02427     }
02428 #ifdef DEBUG
02429   for (unsigned r = 0; r<h; r++)
02430     for (unsigned c =0; c<w; c++)
02431     {
02432       vcl_complex<double> res = fft_matrix[r][c];
02433       vcl_cout << res << '\n';
02434     }
02435 #endif
02436 
02437   brip_vil_float_ops::fft_2d(fft_matrix, w, h, 1);
02438   brip_vil_float_ops::ftt_fourier_2d_reorder(fft_matrix, fourier_matrix);
02439   mag.set_size(w,h);
02440   phase.set_size(w,h);
02441 
02442   // extract magnitude and phase
02443   for (unsigned r = 0; r<h; r++)
02444     for (unsigned c = 0; c<w; c++)
02445     {
02446       float re = (float)fourier_matrix[r][c].real(),
02447         im = (float)fourier_matrix[r][c].imag();
02448       mag(c,r) = vcl_sqrt(re*re + im*im);
02449       phase(c,r) = vcl_atan2(im, re);
02450     }
02451   return true;
02452 }
02453 
02454 bool brip_vil_float_ops::
02455 inverse_fourier_transform(vil_image_view<float> const& mag,
02456                           vil_image_view<float> const& phase,
02457                           vil_image_view<float>& output)
02458 {
02459   unsigned w = mag.ni(), h = mag.nj();
02460   vnl_matrix<vcl_complex<double> > fft_matrix(h, w), fourier_matrix(h, w);
02461   for (unsigned y = 0; y<h; y++)
02462     for (unsigned x =0; x<w; x++)
02463     {
02464       float m = mag(x,y);
02465       float p = phase(x,y);
02466       vcl_complex<double> cv(m*vcl_cos(p), m*vcl_sin(p));
02467       fourier_matrix.put(y, x, cv);
02468     }
02469 
02470   brip_vil_float_ops::ftt_fourier_2d_reorder(fourier_matrix, fft_matrix);
02471   brip_vil_float_ops::fft_2d(fft_matrix, w, h, -1);
02472 
02473   output.set_size(w,h);
02474 
02475   for (unsigned y = 0; y<h; y++)
02476     for (unsigned x = 0; x<w; x++)
02477       output(x,y) = (float)fft_matrix[y][x].real();
02478   return true;
02479 }
02480 
02481 void brip_vil_float_ops::resize(vil_image_view<float> const& input,
02482                                 unsigned width, unsigned height,
02483                                 vil_image_view<float>& output)
02484 {
02485   unsigned w = input.ni(), h = input.nj();
02486   output.set_size(width, height);
02487   for (unsigned y = 0; y<height; y++)
02488     for (unsigned x = 0; x<width; x++)
02489       if (x<w && y<h)
02490         output(x,y) = input(x,y);
02491       else
02492         output(x,y) = 0; // pad with zeroes
02493 }
02494 
02495 //: resize the input to the closest power of two image dimensions
02496 bool brip_vil_float_ops::
02497 resize_to_power_of_two(vil_image_view<float> const& input,
02498                        vil_image_view<float>& output)
02499 {
02500   unsigned max_exp = 13; // we wouldn't want to have such large images in memory
02501   unsigned w = input.ni(), h = input.nj();
02502   unsigned prodw = 1, prodh = 1;
02503   // Find power of two width
02504   unsigned nw, nh;
02505   for (nw = 1; nw<=max_exp; nw++)
02506     if (prodw>w)
02507       break;
02508     else
02509       prodw *= 2;
02510   if (nw==max_exp)
02511     return false;
02512   // Find power of two height
02513   for (nh = 1; nh<=max_exp; nh++)
02514     if (prodh>h)
02515       break;
02516     else
02517       prodh *= 2;
02518   if (nh==max_exp)
02519     return false;
02520   brip_vil_float_ops::resize(input, prodw, prodh, output);
02521 
02522   return true;
02523 }
02524 
02525 //
02526 //: block a periodic signal by suppressing two Gaussian lobes in the frequency domain.
02527 //  The lobes are on the line defined by dir_fx and dir_fy through the
02528 //  dc origin, assumed (0, 0).  The center frequency, f0, is the distance along
02529 //  the line to the center of each blocking lobe (+- f0). radius is the
02530 //  standard deviation of each lobe. Later we can define a "filter" class.
02531 //
02532 float brip_vil_float_ops::gaussian_blocking_filter(float dir_fx,
02533                                                    float dir_fy,
02534                                                    float f0,
02535                                                    float radius,
02536                                                    float fx,
02537                                                    float fy)
02538 {
02539   // normalize dir_fx and dir_fy
02540   float mag = vcl_sqrt(dir_fx*dir_fx + dir_fy*dir_fy);
02541   if (!mag)
02542     return 0.f;
02543   float r2 = 2.f*radius*radius;
02544   float dx = dir_fx/mag, dy = dir_fy/mag;
02545   // compute the centers of each lobe
02546   float fx0p = dx*f0, fy0p = dy*f0;
02547   float fx0m = -dx*f0, fy0m = -dy*f0;
02548   // the squared distance of fx, fy from each center
02549   float d2p = (fx-fx0p)*(fx-fx0p) + (fy-fy0p)*(fy-fy0p);
02550   float d2m = (fx-fx0m)*(fx-fx0m) + (fy-fy0m)*(fy-fy0m);
02551   // use the closest lobe
02552   float d = d2p;
02553   if (d2m<d2p)
02554     d = d2m;
02555   // the gaussian blocking function
02556   float gb = 1.f-(float)vcl_exp(-d/r2);
02557   return gb;
02558 }
02559 
02560 bool brip_vil_float_ops::
02561 spatial_frequency_filter(vil_image_view<float> const& input,
02562                          float dir_fx, float dir_fy,
02563                          float f0, float radius,
02564                          bool output_fourier_mag,
02565                          vil_image_view<float>& output)
02566 {
02567   // Compute the fourier transform of the image.
02568   vil_image_view<float> pow_two, mag, bmag, phase, pow_two_filt;
02569   brip_vil_float_ops::resize_to_power_of_two(input, pow_two);
02570   unsigned Nfx = pow_two.ni(), Nfy = pow_two.nj();
02571 
02572   if (!brip_vil_float_ops::fourier_transform(pow_two, mag, phase))
02573     return false;
02574   bmag.set_size(Nfx, Nfy);
02575 
02576   // filter the magnitude function
02577   float Ofx = Nfx*0.5f, Ofy = Nfy*0.5f;
02578   for (unsigned fy =0; fy<Nfy; fy++)
02579     for (unsigned fx =0; fx<Nfx; fx++)
02580     {
02581       float gb = gaussian_blocking_filter(dir_fx, dir_fy, f0,
02582                                           radius,
02583                                           (float)fx-Ofx, (float)fy-Ofy);
02584       bmag(fx,fy) = mag(fx,fy)*gb;
02585     }
02586   if (output_fourier_mag)
02587   {
02588     output = bmag;
02589     return true;
02590   }
02591   // Transform back
02592   pow_two_filt.set_size(Nfx, Nfy);
02593   brip_vil_float_ops::inverse_fourier_transform(bmag, phase, pow_two_filt);
02594 
02595   // Resize to original input size
02596   brip_vil_float_ops::resize(pow_two_filt, input.ni(), input.nj(), output);
02597   return true;
02598 }
02599 
02600 //----------------------------------------------------------------------
02601 //: Bi-linear interpolation on the neighborhood below.
02602 // \verbatim
02603 //      xr
02604 //   yr 0  x
02605 //      x  x
02606 // \endverbatim
02607 //
02608 double brip_vil_float_ops::
02609 bilinear_interpolation(vil_image_view<float> const& input, double x, double y)
02610 {
02611   // check bounds
02612   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
02613   // the pixel containing the interpolated point
02614   int xr = (int)x, yr = (int)y;
02615   double fx = x-xr, fy = y-yr;
02616   if (xr<0||xr>w-2)
02617     return 0.0;
02618   if (yr<0||yr>h-2)
02619     return 0.0;
02620   double int00 = input(xr, yr), int10 = input(xr+1,yr);
02621   double int01 = input(xr, yr+1), int11 = input(xr+1,yr+1);
02622   double int0 = int00 + fy * (int01 - int00);
02623   double int1 = int10 + fy * (int11 - int10);
02624   float val = (float) (int0 + fx * (int1 - int0));
02625   return val;
02626 }
02627 
02628 //: Transform the input to the output by a homography.
02629 //  If the output size is fixed then only the corresponding
02630 //  region of input image space is transformed.
02631 bool brip_vil_float_ops::homography(vil_image_view<float> const& input,
02632                                     vgl_h_matrix_2d<double> const& H,
02633                                     vil_image_view<float>& output,
02634                                     bool output_size_fixed,
02635                                     float output_fill_value)
02636 {
02637   if (!input)
02638     return false;
02639   // smooth the input to condition interpolation
02640   vil_image_view<float> gimage =
02641     brip_vil_float_ops::gaussian(input, 0.5f);
02642 
02643   // First, there is some rather complex bookeeping to insure that
02644   // the input and output image rois are consistent with the homography.
02645 
02646   // the bounding boxes corresponding to input and output rois
02647   // We also construct polygons since homographies turn boxes into arbitrary
02648   // quadrilaterals.
02649   vsol_box_2d_sptr input_roi, output_roi;
02650   vsol_polygon_2d_sptr input_poly, output_poly;
02651   vgl_h_matrix_2d<double> Hinv;
02652   // set up the roi and poly for the input image
02653   unsigned win = gimage.ni(), hin = gimage.nj();
02654   input_roi = new vsol_box_2d();
02655   input_roi->add_point(0, 0);
02656   input_roi->add_point(win, hin);
02657   input_poly = bsol_algs::poly_from_box(input_roi);
02658   // Case I
02659   //  the output image size and input transform can be adjusted
02660   //  to map the transformed image onto the full range
02661   if (!output_size_fixed)
02662   {
02663     if (!bsol_algs::homography(input_poly, H, output_poly))
02664       return false;
02665     vsol_box_2d_sptr temp = output_poly->get_bounding_box();
02666     output.set_size((int)temp->width(), (int)temp->height());
02667     output.fill(output_fill_value);
02668     //offset the transform so the origin is (0,0)
02669     output_roi = new vsol_box_2d();
02670     output_roi->add_point(0, 0);
02671     output_roi->add_point(temp->width(), temp->height());
02672     vgl_h_matrix_2d<double> t; t.set_identity().set_translation(-temp->get_min_x(),-temp->get_min_y());
02673     Hinv = (t*H).get_inverse();
02674   }
02675   else // Case II, the output image size is fixed so we have to find the
02676   {    // inverse mapping of the output roi and intersect with the input roi
02677     // to determine the domain of the mapping
02678     if (!output)
02679       return false;
02680     //The output roi and poly
02681     unsigned wout = output.ni(), hout = output.nj();
02682     output.fill(output_fill_value);
02683     output_roi = new vsol_box_2d();
02684     output_roi->add_point(0, 0);
02685     output_roi->add_point(wout, hout);
02686     output_poly = bsol_algs::poly_from_box(output_roi);
02687 
02688     //Construct the reverse mapping of the output bounds
02689     vsol_polygon_2d_sptr tpoly;
02690     Hinv = H.get_inverse();
02691     if (!bsol_algs::homography(output_poly, Hinv, tpoly))
02692       return false;
02693 
02694     //form the roi corresponding to the inverse mapped output bounds
02695     vsol_box_2d_sptr tbox = tpoly->get_bounding_box();
02696 
02697     //intersect with the input image bounds to get the input roi
02698     vsol_box_2d_sptr temp;
02699     if (!bsol_algs::intersection(tbox, input_roi, temp))
02700       return false;
02701     input_roi = temp;
02702   }
02703   // At this point we have the correct bounds for the input and
02704   // the output image
02705 
02706   // Iterate over the output image space and map the location of each
02707   // pixel into the input image space. Then carry out interpolation to
02708   // get the value of each output pixel
02709 
02710   // Dimensions of the input image
02711   int ailow  = int(input_roi->get_min_x()+0.9999f); // round up to nearest int
02712   int aihigh = int(input_roi->get_max_x());      // round down to nearest int
02713   int ajlow  = int(input_roi->get_min_y()+0.9999f);
02714   int ajhigh = int(input_roi->get_max_y());
02715 
02716   // Dimensions of the output image
02717   int bilow  = int(output_roi->get_min_x()+0.9999f);
02718   int bihigh = int(output_roi->get_max_x());
02719   int bjlow  = int(output_roi->get_min_y()+0.9999f);
02720   int bjhigh = int(output_roi->get_max_y());
02721 
02722   // The inverse transform is used to map backwards from the output
02723   const vnl_matrix_fixed<double,3,3>& Minv = Hinv.get_matrix();
02724 
02725   // Now use Hinv to transform the image
02726   for (int i = bilow; i<bihigh; i++)
02727     for (int j = bjlow; j<bjhigh; j++)
02728     {
02729       // Transform the pixel
02730       float val;
02731       double u = Minv[0][0] * i + Minv[0][1] * j + Minv[0][2];
02732       double v = Minv[1][0] * i + Minv[1][1] * j + Minv[1][2];
02733       double w = Minv[2][0] * i + Minv[2][1] * j + Minv[2][2];
02734       u /= w;
02735       v /= w;
02736 
02737       // Now do linear interpolation
02738       {
02739         int iu = (int) u;
02740         int iv = (int) v;
02741 
02742         if (iu >= ailow && iu < aihigh-1 &&
02743             iv >= ajlow && iv < ajhigh-1)
02744         {
02745           // Get the neighbouring pixels
02746           //      (u  v)    (u+1  v)
02747           //      (u v+1)   (u+1 v+1)
02748           //
02749           double v00 = gimage(iu, iv);
02750           double v01 = gimage(iu, iv+1);
02751           double v10 = gimage(iu+1,iv);
02752           double v11 = gimage(iu+1, iv+1);
02753 
02754           double fu = u - iu;
02755           double fv = v - iv;
02756           double v0 = v00 + fv * (v01 - v00);
02757           double v1 = v10 + fv * (v11 - v10);
02758           val = (float) (v0 + fu * (v1 - v0));
02759           // Set the value
02760           output(i,j) = val;
02761         }
02762       }
02763     }
02764   return true;
02765 }
02766 
02767 //: rotate the input image counter-clockwise about the image origin.
02768 // Demonstrates the use of image homography
02769 vil_image_view<float>
02770 brip_vil_float_ops::rotate(vil_image_view<float> const& input,
02771                            double theta_deg)
02772 {
02773   vil_image_view<float> out;
02774   if (!input)
02775     return out;
02776   double ang = theta_deg;
02777   // map theta_deg to [0 360]
02778   while (ang>360)
02779     ang-=360;
02780   while (ang<0)
02781     ang+=360;
02782   // convert to radians
02783   double deg_to_rad = vnl_math::pi_over_180;
02784   double rang = deg_to_rad*ang;
02785   vgl_h_matrix_2d<double> H;
02786   H.set_identity().set_rotation(rang);
02787   vil_image_view<float> temp;
02788   // The transform is adjusted to map the full input domain onto
02789   // the output image.
02790   if (!brip_vil_float_ops::homography(input, H, temp))
02791     return out;
02792   return temp;
02793 }
02794 
02795 bool brip_vil_float_ops::chip(vil_image_view<float> const& input,
02796                               vsol_box_2d_sptr const& roi,
02797                               vil_image_view<float>& chip)
02798 {
02799   if (!input||!roi)
02800     return false;
02801   int w = static_cast<int>(input.ni()), h = static_cast<int>(input.nj());
02802   int x_min = (int)roi->get_min_x(), y_min = (int)roi->get_min_y();
02803   int x_max = (int)roi->get_max_x(), y_max = (int)roi->get_max_y();
02804   if (x_min<0)
02805     x_min = 0;
02806   if (y_min<0)
02807     y_min = 0;
02808   if (x_max>w-1)
02809     x_max=w-1;
02810   if (y_max>h-1)
02811     y_max=h-1;
02812   int rw = x_max-x_min, rh = y_max-y_min;
02813   if (rw<=0||rh<=0)
02814     return false;
02815   vil_image_view<float> temp(rw+1, rh+1, 1);
02816   for (int y = y_min; y<=y_max; y++)  //!< changed < to <= to include the boundary points too
02817     for (int x =x_min; x<=x_max; x++)
02818       temp(x-x_min, y-y_min) = input(x, y);
02819   chip = temp;
02820   return true;
02821 }
02822 
02823 //: convert image resource to cropped view according to a roi.
02824 bool brip_vil_float_ops::chip(vil_image_resource_sptr const& image,
02825                               brip_roi_sptr const& roi,
02826                               vil_image_resource_sptr & chip)
02827 {
02828   // image bounds
02829   unsigned ni = image->ni(), nj = image->nj();
02830 
02831   // image bounds for the chip
02832   unsigned cm = roi->cmin(0), rm = roi->rmin(0);
02833   unsigned niv = roi->csize(0), njv = roi->rsize(0);
02834   vcl_cout << "org(" << cm << ' ' << rm << "):size(" << niv
02835            << ' ' << njv << ")\n" << vcl_flush;
02836   // check bounds
02837   if (cm>ni-1||rm>nj-1||cm+niv>ni||rm+njv>nj)
02838     return false;
02839   vil_pixel_format pix_format = image->pixel_format();
02840   // get an appropriate image view for scalar images we care about
02841   if (image->nplanes()==1)
02842   {
02843     if (pix_format==VIL_PIXEL_FORMAT_BYTE)
02844     {
02845       vil_image_view<unsigned char> temp = image->get_view(cm, niv, rm, njv);
02846       if (!temp) return false;
02847       chip = vil_new_image_resource_of_view(temp);
02848       return true;
02849     }
02850     else if (pix_format==VIL_PIXEL_FORMAT_UINT_16)
02851     {
02852       vil_image_view<unsigned short> temp = image->get_view(cm, niv, rm, njv);
02853       if (!temp) return false;
02854       chip = vil_new_image_resource_of_view(temp);
02855       vcl_cout << "resc(" << chip->ni() << ' ' << chip->nj()<< ")\n"
02856                << vcl_flush;
02857       return true;
02858     }
02859     else if (pix_format==VIL_PIXEL_FORMAT_FLOAT)
02860     {
02861       vil_image_view<float> temp = image->get_view(cm, niv, rm, njv);
02862       if (!temp) return false;
02863       chip = vil_new_image_resource_of_view(temp);
02864       return true;
02865     }
02866   }
02867 
02868   // color data
02869   if (image->nplanes()==3)
02870   {
02871     if (pix_format==VIL_PIXEL_FORMAT_BYTE) //the only way now
02872     {
02873       //extract view corresponding to region of interest
02874       vil_image_view<vil_rgb<vxl_byte> > temp = image->get_view(cm, niv, rm, njv);
02875       if (!temp) return false;
02876       chip = vil_new_image_resource_of_view(temp);
02877       return true;
02878     }
02879   }
02880   return false;
02881 }
02882 
02883 bool brip_vil_float_ops::
02884 chip(vcl_vector<vil_image_resource_sptr> const& images,
02885      brip_roi_sptr const& roi,
02886      vcl_vector<vil_image_resource_sptr>& chips)
02887 {
02888   for (vcl_vector<vil_image_resource_sptr>::const_iterator iit = images.begin(); iit != images.end(); ++iit) {
02889     vil_image_resource_sptr chip;
02890     if (!brip_vil_float_ops::chip(*iit, roi, chip))
02891       return false;
02892     chips.push_back(chip);
02893   }
02894   return true;
02895 }
02896 
02897 //: perform normalized cross-correlation at a sub-pixel location.
02898 // Thus all the pixel values are interpolated.
02899 float brip_vil_float_ops::
02900 cross_correlate(vil_image_view<float> const& image1,
02901                 vil_image_view<float> const& image2,
02902                 float x, float y, int radius, float intensity_thresh)
02903 {
02904   unsigned w1 = image1.ni(), h1 = image1.nj();
02905   unsigned w2 = image1.ni(), h2 = image1.nj();
02906   // bounds checks
02907   if (w1!=w2||h1!=h2)
02908     return -1;
02909   if (x<radius||x>w1-radius-1||y<radius||y>h1-radius-1)
02910     return -1;
02911 
02912   // accumulate correlation sums,
02913   // bi-linear interpolate the values
02914   int s = 2*radius+1;
02915   double area = s*s;
02916   double sI1=0, sI2=0, sI1I1=0, sI2I2=0, sI1I2=0;
02917   for (int y0 = -10*radius; y0<=10*radius; ++y0)
02918     for (int x0 = -10*radius; x0<=10*radius; ++x0)
02919     {
02920       float xp = x+0.1f*(float)x0, yp = y+0.1f*(float)y0;
02921       double v1 = brip_vil_float_ops::bilinear_interpolation(image1, xp, yp);
02922       double v2 = brip_vil_float_ops::bilinear_interpolation(image2, xp, yp);
02923       sI1 += v1;
02924       sI2 += v2;
02925       sI1I1 += v1*v1;
02926       sI2I2 += v2*v2;
02927       sI1I2 += v1*v2;
02928     }
02929   // compute correlation.
02930   float cc = cross_corr(area, sI1, sI2, sI1I1, sI2I2, sI1I2, intensity_thresh);
02931   return cc;
02932 }
02933 
02934 //: r0 is the image from which to read the new intensity values.
02935 //  \param r is the summing array row in which the values are to be accumulated
02936 static bool update_row(vil_image_view<float> const& image1,
02937                        vil_image_view<float> const& image2,
02938                        int r0, int r,
02939                        vbl_array_2d<double>& SI1,
02940                        vbl_array_2d<double>& SI2,
02941                        vbl_array_2d<double>& SI1I1,
02942                        vbl_array_2d<double>& SI2I2,
02943                        vbl_array_2d<double>& SI1I2)
02944 {
02945   unsigned w1 = image1.ni();
02946   unsigned w2 = image2.ni();
02947   int h1 = image1.nj();
02948   int h2 = image2.nj();
02949   if (w1!=w2||h1!=h2||r<0||r>=h1)
02950     return false;
02951   double i10 = image1(0,r0), i20 = image2(0,r0);
02952   SI1[r][0] = i10; SI2[r][0] = i20; SI1I1[r][0]=i10*i10;
02953   SI2I2[r][0]=i20*i20; SI1I2[r][0]=i10*i20;
02954   for (unsigned c = 1; c<w1; c++)
02955   {
02956     double i1c = image1(c,r0);
02957     double i2c = image2(c,r0);
02958     SI1[r][c]    = SI1[r][c-1]+i1c;
02959     SI2[r][c]    = SI2[r][c-1]+i2c;
02960     SI1I1[r][c]  = SI1I1[r][c-1]+ i1c*i1c;
02961     SI2I2[r][c]  = SI2I2[r][c-1]+ i2c*i2c;
02962     SI1I2[r][c]  = SI1I2[r][c-1]+ i1c*i2c;
02963   }
02964   return true;
02965 }
02966 
02967 static bool initialize_slice(vil_image_view<float> const& image1,
02968                              vil_image_view<float> const& image2,
02969                              unsigned radius,
02970                              vbl_array_2d<double>& SI1,
02971                              vbl_array_2d<double>& SI2,
02972                              vbl_array_2d<double>& SI1I1,
02973                              vbl_array_2d<double>& SI2I2,
02974                              vbl_array_2d<double>& SI1I2)
02975 {
02976   for (unsigned r = 0; r<=2*radius; r++)
02977     if (!update_row(image1, image2, r, r, SI1, SI2, SI1I1, SI2I2, SI1I2))
02978       return false;
02979   return true;
02980 }
02981 
02982 static bool collapse_slice(vbl_array_2d<double> const& SI1,
02983                            vbl_array_2d<double> const& SI2,
02984                            vbl_array_2d<double> const& SI1I1,
02985                            vbl_array_2d<double> const& SI2I2,
02986                            vbl_array_2d<double> const& SI1I2,
02987                            vbl_array_1d<double>& dSI1,
02988                            vbl_array_1d<double>& dSI2,
02989                            vbl_array_1d<double>& dSI1I1,
02990                            vbl_array_1d<double>& dSI2I2,
02991                            vbl_array_1d<double>& dSI1I2)
02992 {
02993   // sanity check
02994   unsigned w = SI1.cols(), h = SI1.rows();
02995   unsigned dw = dSI1.size();
02996   if (dw!=w)
02997     return false;
02998 
02999   for (unsigned c = 0; c<w; c++)
03000   {
03001     dSI1[c]=0; dSI2[c]=0; dSI1I1[c]=0;
03002     dSI2I2[c]=0; dSI1I2[c]=0;
03003     for (unsigned r = 0; r<h; r++)
03004     {
03005       dSI1[c] += SI1[r][c];
03006       dSI2[c] += SI2[r][c];
03007       dSI1I1[c] += SI1I1[r][c];
03008       dSI2I2[c] += SI2I2[r][c];
03009       dSI1I2[c] += SI1I2[r][c];
03010     }
03011   }
03012   return true;
03013 }
03014 
03015 static bool cross_correlate_row(int radius,
03016                                 vbl_array_1d<double>& dSI1,
03017                                 vbl_array_1d<double>& dSI2,
03018                                 vbl_array_1d<double>& dSI1I1,
03019                                 vbl_array_1d<double>& dSI2I2,
03020                                 vbl_array_1d<double>& dSI1I2,
03021                                 float intensity_thresh,
03022                                 vbl_array_1d<float>& cc)
03023 {
03024   // sanity check
03025   unsigned w = dSI1.size(), wc = cc.size();
03026   if (!w||!wc||w!=wc)
03027     return false;
03028   int s = 2*radius+1;
03029   double area = s*s;
03030   // the general case
03031   double si1=dSI1[s-1], si2=dSI2[s-1], si1i1=dSI1I1[s-1],
03032     si2i2=dSI2I2[s-1], si1i2=dSI1I2[s-1];
03033   cc[radius]= cross_corr(area, si1, si2, si1i1, si2i2, si1i2, intensity_thresh);
03034   // the remaining columns
03035   for (unsigned c = radius+1; c+radius<w; c++)
03036   {
03037     si1=dSI1[c+radius]-dSI1[c-radius-1];
03038     si2=dSI2[c+radius]-dSI2[c-radius-1];
03039     si1i1=dSI1I1[c+radius]-dSI1I1[c-radius-1];
03040     si2i2=dSI2I2[c+radius]-dSI2I2[c-radius-1];
03041     si1i2=dSI1I2[c+radius]-dSI1I2[c-radius-1];
03042     cc[c] = cross_corr(area, si1, si2, si1i1, si2i2, si1i2, intensity_thresh);
03043   }
03044   return true;
03045 }
03046 
03047 static void advance_rows(vbl_array_2d<double>& S)
03048 {
03049   unsigned nr = S.rows(), nc = S.cols();
03050   for (unsigned r = 0; r<nr-1; r++)
03051     for (unsigned c =0; c<nc; c++)
03052       S[r][c]=S[r+1][c];
03053 }
03054 
03055 static bool output_cc_row(int r0, vbl_array_1d<float> const& cc,
03056                           vil_image_view<float>& out)
03057 {
03058   unsigned n = cc.size(), w = out.ni();
03059   if (n!=w)
03060     return false;
03061   for (unsigned c = 0; c<w; c++)
03062     out(c, r0) = cc[c];
03063   return true;
03064 }
03065 
03066 bool brip_vil_float_ops::
03067 cross_correlate(vil_image_view<float> const& image1,
03068                 vil_image_view<float> const& image2,
03069                 vil_image_view<float>& out,
03070                 int radius, float intensity_thresh)
03071 {
03072   vul_timer t;
03073   unsigned w = image1.ni(), h = image1.nj();
03074   unsigned w2 = image2.ni(), h2 = image2.nj();
03075   // sizes must match
03076   if (w!=w2||h!=h2)
03077   {
03078     vcl_cout << "In brip_vil_float_ops::cross_correlate(..) -"
03079              << " image sizes don't match\n";
03080     return out;
03081   }
03082   out.set_size(w, h);
03083   out.fill(0.f);
03084   int s = 2*radius+1;
03085   // Create the running sum slices
03086   vbl_array_2d<double> SI1(s,w), SI2(s,w),
03087     SI1I1(s,w), SI2I2(s,w), SI1I2(s,w);
03088   vbl_array_1d<float> cc(w, 0.f);
03089   vbl_array_1d<double> dSI1(w, 0.0), dSI2(w, 0.0),
03090     dSI1I1(w, 0.0), dSI2I2(w, 0.0), dSI1I2(w, 0.0);
03091   initialize_slice(image1, image2, radius, SI1, SI2, SI1I1, SI2I2, SI1I2);
03092   if (!collapse_slice(SI1,  SI2,  SI1I1,  SI2I2,  SI1I2,
03093                       dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2))
03094     return false;
03095   unsigned r0 = radius;
03096   for (; r0+radius+1<h; r0++)
03097   {
03098     if (r0==5)
03099       r0=5;
03100 #ifdef DEBUG
03101     vcl_cout << "r0 " << r0 << '\n';
03102 #endif
03103     if (!cross_correlate_row(radius, dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2,
03104                              intensity_thresh, cc))
03105       return false;
03106 #ifdef DEBUG
03107     vcl_cout << '\n';
03108 #endif
03109     advance_rows(SI1); advance_rows(SI2);  advance_rows(SI1I1);
03110     advance_rows(SI2I2); advance_rows(SI1I2);
03111     if (!update_row(image1, image2, r0+radius+1, 2*radius,
03112                     SI1, SI2, SI1I1, SI2I2, SI1I2))
03113       return false;
03114     if (!collapse_slice(SI1,  SI2,  SI1I1,  SI2I2,  SI1I2,
03115                         dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2))
03116       return false;
03117     if (!output_cc_row(r0, cc, out))
03118       return out;
03119   }
03120   // handle the last row
03121 #ifdef DEBUG
03122   vcl_cout << "r0 " << r0 << '\n';
03123 #endif
03124   if (!cross_correlate_row(radius, dSI1, dSI2, dSI1I1, dSI2I2, dSI1I2,
03125                            intensity_thresh, cc))
03126     return false;
03127 #ifdef DEBUG
03128   vcl_cout << '\n';
03129 #endif
03130   if (!output_cc_row(r0, cc, out))
03131     return false;
03132 #ifdef DEBUG
03133   vcl_cout << "RunningSumCrossCorrelation for " << w*h/1000.0f << " k pixels in "
03134            << t.real() << " msecs\n"<< vcl_flush;
03135 #endif
03136   return true;
03137 }
03138 
03139 vil_image_resource_sptr
03140 brip_vil_float_ops::sum(vil_image_resource_sptr const& img0,
03141                         vil_image_resource_sptr const& img1)
03142 {
03143   vil_image_view<float> op0 = brip_vil_float_ops::convert_to_float(img0);
03144   vil_image_view<float> op1 = brip_vil_float_ops::convert_to_float(img1);
03145   vil_image_view<float> sum;
03146   vil_math_image_sum(op0, op1, sum);
03147 
03148   // find out the types of the input images for now, only do greyscale operands
03149   vil_pixel_format pix_format0 = img0->pixel_format();
03150   vil_pixel_format pix_format1 = img1->pixel_format();
03151   if (pix_format0 == VIL_PIXEL_FORMAT_FLOAT ||
03152       pix_format1 == VIL_PIXEL_FORMAT_FLOAT)
03153     {
03154       return vil_new_image_resource_of_view(sum);
03155     }
03156 
03157   if (pix_format0 == VIL_PIXEL_FORMAT_UINT_16 ||
03158       pix_format1 == VIL_PIXEL_FORMAT_UINT_16)
03159     {
03160       vil_image_view<vxl_uint_16> res =
03161         brip_vil_float_ops::convert_to_short(vil_new_image_resource_of_view(sum));
03162       return vil_new_image_resource_of_view(res);
03163     }
03164 
03165   if (pix_format0 == VIL_PIXEL_FORMAT_BYTE ||
03166       pix_format1 == VIL_PIXEL_FORMAT_BYTE)
03167     {
03168       vil_image_view<vxl_byte> res =
03169         brip_vil_float_ops::convert_to_byte(sum);
03170       return vil_new_image_resource_of_view(res);
03171     }
03172   // for color return the float image
03173   return vil_new_image_resource_of_view(sum);
03174 }
03175 
03176 // Compute img0 - img1
03177 vil_image_resource_sptr
03178 brip_vil_float_ops::difference(vil_image_resource_sptr const& img0,
03179                                vil_image_resource_sptr const& img1)
03180 {
03181   vil_image_view<float> op0 = brip_vil_float_ops::convert_to_float(img0);
03182   vil_image_view<float> op1 = brip_vil_float_ops::convert_to_float(img1);
03183   vil_image_view<float> diff;
03184   vil_math_image_difference(op0, op1, diff);
03185 
03186   // find out the types of the input images for now, only do greyscale operands
03187   vil_pixel_format pix_format0 = img0->pixel_format();
03188   vil_pixel_format pix_format1 = img1->pixel_format();
03189   if (pix_format0 == VIL_PIXEL_FORMAT_FLOAT ||
03190       pix_format1 == VIL_PIXEL_FORMAT_FLOAT)
03191     {
03192       return vil_new_image_resource_of_view(diff);
03193     }
03194 
03195   if (pix_format0 == VIL_PIXEL_FORMAT_UINT_16 ||
03196       pix_format1 == VIL_PIXEL_FORMAT_UINT_16)
03197     {
03198       vil_image_view<vxl_uint_16> res =
03199         brip_vil_float_ops::convert_to_short(vil_new_image_resource_of_view(diff));
03200       return vil_new_image_resource_of_view(res);
03201     }
03202 
03203   if (pix_format0 == VIL_PIXEL_FORMAT_BYTE ||
03204       pix_format1 == VIL_PIXEL_FORMAT_BYTE)
03205     {
03206       vil_image_view<vxl_byte> res = brip_vil_float_ops::convert_to_byte(diff);
03207       return vil_new_image_resource_of_view(res);
03208     }
03209   // for color return the float image
03210   return vil_new_image_resource_of_view(diff);
03211 }
03212 
03213 //: Compute the entropy of the intensity of a region
03214 //  Note no bounds checking!
03215 float brip_vil_float_ops::entropy_i(const unsigned i, const unsigned j,
03216                                     const unsigned i_radius,
03217                                     const unsigned j_radius,
03218                                     vil_image_view<float> const& intensity,
03219                                     const float range, const unsigned bins)
03220 {
03221   bsta_histogram<float> hi(range, bins);
03222   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03223   for (int dj = -jr; dj<=jr; ++dj)
03224     for (int di = -ir; di<=ir; ++di)
03225     {
03226       float inten = intensity(i+di, j+dj);
03227       hi.upcount(inten, 1.0f);
03228     }
03229   return hi.entropy();
03230 }
03231 
03232 //: Compute the entropy of the gradient direction of a region
03233 //  Note no bounds checking!
03234 float brip_vil_float_ops::entropy_g(const unsigned i, const unsigned j,
03235                                     const unsigned i_radius,
03236                                     const unsigned j_radius,
03237                                     vil_image_view<float> const& gradx,
03238                                     vil_image_view<float> const& grady,
03239                                     const float range, const unsigned bins)
03240 {
03241   bsta_histogram<float> hg(range, bins);
03242   static const float deg_rad = (float)(vnl_math::deg_per_rad);
03243   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03244   for (int dj = -jr; dj<=jr; ++dj)
03245     for (int di = -ir; di<=ir; ++di)
03246     {
03247       float Ix = gradx(i+di, j+dj), Iy = grady(i+di, j+dj);
03248       float ang = deg_rad*vcl_atan2(Iy, Ix) + 180.0f;
03249       float mag = vcl_abs(Ix)+vcl_abs(Iy);
03250       hg.upcount(ang, mag);
03251     }
03252   return hg.entropy();
03253 }
03254 
03255 //: Compute the hue and saturation entropy of a region about the specified pixel
03256 //  Note no bounds checking!
03257 float brip_vil_float_ops::entropy_hs(const unsigned i, const unsigned j,
03258                                      const unsigned i_radius,
03259                                      const unsigned j_radius,
03260                                      vil_image_view<float> const& hue,
03261                                      vil_image_view<float> const& sat,
03262                                      const float range, const unsigned bins)
03263 {
03264   bsta_histogram<float> hg(range, bins);
03265   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03266   for (int dj = -jr; dj<=jr; ++dj)
03267     for (int di = -ir; di<=ir; ++di)
03268     {
03269       float h = hue(i+di, j+dj), s = sat(i+di, j+dj);
03270       hg.upcount(h, s);
03271     }
03272   return hg.entropy();
03273 }
03274 
03275 vil_image_view<float>
03276 brip_vil_float_ops::entropy(const unsigned i_radius,
03277                             const unsigned j_radius,
03278                             const unsigned step,
03279                             vil_image_resource_sptr const& img,
03280                             const float sigma,
03281                             const unsigned bins,
03282                             const bool intensity,
03283                             const bool gradient,
03284                             const bool ihs)
03285 {
03286   vil_image_view<float> ent;
03287   if (!intensity&&!gradient&&!ihs)
03288   {
03289     vcl_cout << "In brip_vil_float_ops::entropy(.) - No computation to do\n";
03290     return ent;
03291   }
03292 
03293   vil_image_view<float> fimage = brip_vil_float_ops::convert_to_float(img);
03294   vil_image_view<float> gimage =
03295     brip_vil_float_ops::gaussian(fimage, sigma);
03296 
03297   unsigned ni = img->ni(), nj = img->nj();
03298   ent.set_size(ni/step+1, nj/step+1);
03299   ent.fill(0.0f);
03300   if (intensity)
03301     for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
03302       for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
03303         ent(i/step,j/step) =
03304           brip_vil_float_ops::entropy_i(i, j, i_radius, j_radius, gimage, 255.0f, bins);
03305 
03306   if (gradient)
03307   {
03308     vil_image_view<float> grad_x, grad_y;
03309     grad_x.set_size(ni, nj);
03310     grad_y.set_size(ni, nj);
03311     brip_vil_float_ops::gradient_3x3 (gimage , grad_x , grad_y);
03312     for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
03313       for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
03314         ent(i/step,j/step) +=
03315           brip_vil_float_ops::entropy_g(i, j, i_radius, j_radius,
03316                                         grad_x, grad_y);
03317   }
03318   if (ihs&&img->nplanes()==3)
03319   {
03320     vil_image_view<float> inten, hue, sat;
03321     vil_image_view<vil_rgb<vxl_byte> > cimage = img->get_view();
03322     brip_vil_float_ops::convert_to_IHS(cimage, inten, hue, sat);
03323     for (unsigned j = j_radius; j<(nj-j_radius); j+=step)
03324       for (unsigned i = i_radius; i<(ni-i_radius); i+=step)
03325         ent(i/step,j/step) +=
03326           brip_vil_float_ops::entropy_hs(i, j, i_radius, j_radius,
03327                                          hue, sat);
03328   }
03329   return ent;
03330 }
03331 
03332 float brip_vil_float_ops::minfo_i(const unsigned i0, const unsigned j0,
03333                                   const unsigned i1, const unsigned j1,
03334                                   const unsigned i_radius,
03335                                   const unsigned j_radius,
03336                                   vil_image_view<float> const& int0,
03337                                   vil_image_view<float> const& int1,
03338                                   const float range, const unsigned bins)
03339 {
03340   bsta_histogram<float> hi0(range, bins);
03341   bsta_histogram<float> hi1(range, bins);
03342   bsta_joint_histogram<float> hji(range, bins);
03343   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03344   for (int dj = -jr; dj<=jr; ++dj)
03345     for (int di = -ir; di<=ir; ++di)
03346     {
03347       float inten0 = int0(i0+di, j0+dj);
03348       float inten1 = int1(i1+di, j1+dj);
03349       hi0.upcount(inten0, 1.0f);
03350       hi1.upcount(inten1, 1.0f);
03351       hji.upcount(inten0, 1.0f, inten1, 1.0f);
03352     }
03353   float H0 = hi0.entropy();
03354   float H1 = hi1.entropy();
03355   float HJ = hji.entropy();
03356   float minfo_i = H0 + H1 - HJ;
03357 #ifdef DEBUG
03358   if (minfo<0)
03359     vcl_cout << "intensity MI LT 0 " << minfo << vcl_endl;
03360 #endif
03361   return minfo_i;
03362 }
03363 
03364 float brip_vil_float_ops::minfo_g(const unsigned i0, const unsigned j0,
03365                                   const unsigned i1, const unsigned j1,
03366                                   const unsigned i_radius,
03367                                   const unsigned j_radius,
03368                                   vil_image_view<float> const& gradx0,
03369                                   vil_image_view<float> const& grady0,
03370                                   vil_image_view<float> const& gradx1,
03371                                   vil_image_view<float> const& grady1,
03372                                   const float range, const unsigned bins)
03373 {
03374   bsta_histogram<float> hg0(range, bins);
03375   bsta_histogram<float> hg1(range, bins);
03376   bsta_joint_histogram<float> hjg(range, bins);
03377   static const float deg_rad = (float)(vnl_math::deg_per_rad);
03378   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03379   for (int dj = -jr; dj<=jr; ++dj)
03380     for (int di = -ir; di<=ir; ++di)
03381     {
03382       float Ix0 = gradx0(i0+di, j0+dj), Iy0 = grady0(i0+di, j0+dj);
03383       float Ix1 = gradx1(i1+di, j1+dj), Iy1 = grady1(i1+di, j1+dj);
03384       float ang0 = deg_rad*vcl_atan2(Iy0, Ix0) + 180.0f;
03385       float ang1 = deg_rad*vcl_atan2(Iy1, Ix1) + 180.0f;
03386       float mag0 = vcl_abs(Ix0)+vcl_abs(Iy0);
03387       float mag1 = vcl_abs(Ix1)+vcl_abs(Iy1);
03388       hg0.upcount(ang0, mag0);
03389       hg1.upcount(ang1, mag1);
03390       hjg.upcount(ang0, mag0, ang1, mag1);
03391     }
03392   float H0 = hg0.entropy();
03393   float H1 = hg1.entropy();
03394   float HJ = hjg.entropy();
03395   float minfo_g = H0 + H1 - HJ;
03396 #ifdef DEBUG
03397   if (minfo<0)
03398     vcl_cout << "gradient MI LT 0 " << minfo << vcl_endl;
03399 #endif
03400   return minfo_g;
03401 }
03402 
03403 float brip_vil_float_ops::minfo_hs(const unsigned i0, const unsigned j0,
03404                                    const unsigned i1, const unsigned j1,
03405                                    const unsigned i_radius,
03406                                    const unsigned j_radius,
03407                                    vil_image_view<float> const& hue0,
03408                                    vil_image_view<float> const& sat0,
03409                                    vil_image_view<float> const& hue1,
03410                                    vil_image_view<float> const& sat1,
03411                                    const float range, const unsigned bins)
03412 {
03413   bsta_histogram<float> hh0(range, bins);
03414   bsta_histogram<float> hh1(range, bins);
03415   bsta_joint_histogram<float> hjh(range, bins);
03416   int ir = static_cast<int>(i_radius), jr = static_cast<int>(j_radius);
03417   for (int dj = -jr; dj<=jr; ++dj)
03418     for (int di = -ir; di<=ir; ++di)
03419     {
03420       float h0 = hue0(i0+di, j0+dj), s0 = sat0(i0+di, j0+dj);
03421       float h1 = hue1(i1+di, j1+dj), s1 = sat1(i1+di, j1+dj);
03422       hh0.upcount(h0, s0);
03423       hh1.upcount(h1, s1);
03424       hjh.upcount(h0, s0, h1, s1);
03425     }
03426   float H0 = hh0.entropy();
03427   float H1 = hh1.entropy();
03428   float HJ = hjh.entropy();
03429   float minfo_h = H0 + H1 - HJ;
03430 #ifdef DEBUG
03431   if (minfo<0)
03432     vcl_cout << "color MI LT 0 " << minfo << vcl_endl;
03433 #endif
03434   return minfo_h;
03435 }
03436 
03437 bool brip_vil_float_ops::minfo(const unsigned i_radius,
03438                                const unsigned j_radius,
03439                                const unsigned step,
03440                                vil_image_resource_sptr const& img0,
03441                                vil_image_resource_sptr const& img1,
03442                                vil_image_view<float>& MI0,
03443                                vil_image_view<float>& MI1,
03444                                const float sigma,
03445                                const bool intensity,
03446                                const bool gradient,
03447                                const bool ihs)
03448 {
03449   if (!intensity&&!gradient&&!ihs)
03450   {
03451     vcl_cout << "In brip_vil_float_ops::minforopy(.) - No computation to do\n";
03452     return false;
03453   }
03454 
03455   vil_image_view<float> fimage0 = brip_vil_float_ops::convert_to_float(img0);
03456   vil_image_view<float> gimage0 =
03457     brip_vil_float_ops::gaussian(fimage0, sigma);
03458 
03459   vil_image_view<float> fimage1 = brip_vil_float_ops::convert_to_float(img1);
03460   vil_image_view<float> gimage1 =
03461     brip_vil_float_ops::gaussian(fimage1, sigma);
03462 
03463   unsigned ni0 = img0->ni(), nj0 = img0->nj();
03464   unsigned ni1 = img1->ni(), nj1 = img1->nj();
03465   unsigned ilimit = 2*i_radius +1, jlimit = 2*j_radius+1;
03466   if (ni0<ilimit||nj0<jlimit||ni1<ilimit||nj1<jlimit)
03467   {
03468     vcl_cout << "In brip_vil_float_ops::minfo(...) - image too small\n";
03469     return false;
03470   }
03471   MI0.set_size(ni0/step+1, nj0/step+1); MI0.fill(0.0f);
03472   MI1.set_size(ni1/step+1, nj1/step+1); MI1.fill(0.0f);
03473   if (intensity)
03474     for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
03475       for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
03476         for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
03477           for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
03478           {
03479             float minfo = brip_vil_float_ops::minfo_i(i0, j0,i1, j1,
03480                                                       i_radius, j_radius,
03481                                                       gimage0, gimage1);
03482             MI0(i0/step,j0/step) = minfo;
03483             MI1(i1/step,j1/step) = minfo;
03484           }
03485   if (gradient)
03486   {
03487     vil_image_view<float> grad_x0, grad_y0, grad_x1, grad_y1;
03488     grad_x0.set_size(ni0, nj0);
03489     grad_y0.set_size(ni0, nj0);
03490     grad_x1.set_size(ni1, nj1);
03491     grad_y1.set_size(ni1, nj1);
03492     brip_vil_float_ops::gradient_3x3 (gimage0 , grad_x0 , grad_y0);
03493     brip_vil_float_ops::gradient_3x3 (gimage1 , grad_x1 , grad_y1);
03494     for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
03495       for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
03496         for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
03497           for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
03498           {
03499             float minfo = brip_vil_float_ops::minfo_g(i0, j0,i1, j1,
03500                                                       i_radius, j_radius,
03501                                                       grad_x0, grad_y0,
03502                                                       grad_x1, grad_y1);
03503             MI0(i0/step,j0/step) += minfo;
03504             MI1(i1/step,j1/step) += minfo;
03505           }
03506   }
03507   if (ihs&&img0->nplanes()==3&&img1->nplanes()==3)
03508   {
03509     vil_image_view<float> inten0, hue0, sat0;
03510     vil_image_view<float> inten1, hue1, sat1;
03511     vil_image_view<vil_rgb<vxl_byte> > cimage0 = img0->get_view();
03512     vil_image_view<vil_rgb<vxl_byte> > cimage1 = img1->get_view();
03513     brip_vil_float_ops::convert_to_IHS(cimage0, inten0, hue0, sat0);
03514     brip_vil_float_ops::convert_to_IHS(cimage1, inten1, hue1, sat1);
03515 
03516     for (unsigned j0 = j_radius; j0<(nj0-j_radius); j0+=step)
03517       for (unsigned i0 = i_radius; i0<(ni0-i_radius); i0+=step)
03518         for (unsigned j1 = j_radius; j1<(nj1-j_radius); j1+=step)
03519           for (unsigned i1 = i_radius; i1<(ni1-i_radius); i1+=step)
03520           {
03521             float minfo = brip_vil_float_ops::minfo_hs(i0, j0,i1, j1,
03522                                                        i_radius, j_radius,
03523                                                        hue0, sat0,
03524                                                        hue1, sat1);
03525             MI0(i0/step,j0/step) += minfo;
03526             MI1(i1/step,j1/step) += minfo;
03527           }
03528   }
03529   return true;
03530 }
03531 
03532 // compute the average of the image intensity within the specified region
03533 float brip_vil_float_ops::
03534 average_in_box(vil_image_view<float> const& v, vgl_box_2d<double> const&  box)
03535 {
03536   vgl_point_2d<double> p0 = box.min_point();
03537   unsigned i0 = static_cast<unsigned>(p0.x()), j0 = static_cast<unsigned>(p0.y());
03538   vgl_point_2d<double> p1 = box.max_point();
03539   unsigned i1 = static_cast<unsigned>(p1.x()), j1 = static_cast<unsigned>(p1.y());
03540   float n = 0.0f;
03541   float sum = 0.0f;
03542   for (unsigned i = i0; i<=i1; ++i)
03543     for (unsigned j = j0; j<=j1; ++j, ++n)
03544       sum += v(i,j);
03545   if (n>0)
03546     sum /= n;
03547   return sum;
03548 }
03549 
03550 #if 0 // For now remove dependency on vimt. Save for illustration
03551 bool brip_vil_float_ops::vimt_homography(vil_image_view<float> const& curr_view,
03552                                          vgl_h_matrix_2d<double>const& H,
03553                                          vil_image_view<float>& output)
03554 {
03555   int bimg_ni;
03556   int bimg_nj;
03557 
03558   int offset_i;
03559   int offset_j;
03560 
03561   vbl_bounding_box<double,2> box;
03562 
03563   unsigned ni =  curr_view.ni(), nj =  curr_view.nj();
03564   vimt_transform_2d p;
03565   vnl_double_3x3 Mh = H.get_matrix();
03566   vnl_double_2x3 M;
03567   for (unsigned c = 0; c<3; ++c)
03568     for (unsigned r = 0; r<2; ++r)
03569       M[r][c] = Mh[r][c]/Mh[2][2];
03570   p.set_affine(M);
03571   box.update(p(0,0).x(),p(0,0).y());
03572   box.update(p(0,nj).x(),p(0,nj).y());
03573   box.update(p(ni,0).x(),p(ni,0).y());
03574   box.update(p(ni,nj).x(),p(ni,nj).y());
03575 
03576   bimg_ni=(int)vcl_ceil(box.max()[0]-box.min()[0]);
03577   bimg_nj=(int)vcl_ceil(box.max()[1]-box.min()[1]);
03578 
03579   offset_i=(int)vcl_ceil(0-box.min()[0]);
03580   offset_j=(int)vcl_ceil(0-box.min()[1]);
03581 
03582   vimt_image_2d_of<float> sample_im;
03583   vgl_point_2d<double> q(-offset_i,-offset_j);
03584   vgl_vector_2d<double> u(1,0);
03585   vgl_vector_2d<double> v(0,1);
03586 
03587   vimt_image_2d_of<float> curr_img(curr_view,p);
03588   vimt_resample_bilin(curr_img,sample_im,q,u,v,bimg_ni,bimg_nj);
03589   output = sample_im.image();
03590   return true;
03591 }
03592 
03593 #endif // 0
03594 
03595 vcl_vector<float> brip_vil_float_ops::scan_region(vil_image_resource_sptr img,
03596                                                   vgl_polygon<double> poly,
03597                                                   float& min,
03598                                                   float& max)
03599 {
03600   vcl_vector<float> pixels;
03601   min = vnl_numeric_traits<float>::maxval;
03602   max = 0;
03603   unsigned np = img->nplanes();
03604   vgl_polygon_scan_iterator<double> si(poly, false);
03605   if (img->pixel_format()==VIL_PIXEL_FORMAT_BYTE)
03606   {
03607     if (np==1) // single plane
03608     {
03609       for (si.reset(); si.next();)
03610       {
03611         unsigned j = static_cast<unsigned>(si.scany());
03612         for (int x = si.startx(); x<=si.endx(); ++x)
03613         {
03614           unsigned i = static_cast<unsigned>(x);
03615           vil_image_view<unsigned char> v = img->get_view(i, 1,j, 1);
03616           float fv = static_cast<float>(v(0,0));
03617           if (fv<min) min = fv;
03618           if (fv>max) max = fv;
03619           pixels.push_back(fv);
03620         }
03621       }
03622       return pixels;
03623     }
03624     else
03625     {
03626       for (si.reset(); si.next();)
03627       {
03628         unsigned j = static_cast<unsigned>(si.scany());
03629         for (int x = si.startx(); x<=si.endx(); ++x)
03630         {
03631           unsigned i = static_cast<unsigned>(x);
03632           vil_image_view<unsigned char> v = img->get_view(i, 1,j, 1);
03633           float fv = 0;
03634           for (unsigned p = 0; p<np; ++p)
03635             fv += v(0,0,p);
03636           fv/=3;
03637           if (fv<min) min = fv;
03638           if (fv>max) max = fv;
03639           pixels.push_back(fv);
03640         }
03641       }
03642     }
03643     return pixels;
03644   }
03645   else if (img->pixel_format()==VIL_PIXEL_FORMAT_UINT_16)
03646   {
03647     if (np) // single plane
03648     {
03649       for (si.reset(); si.next();)
03650       {
03651         unsigned j = static_cast<unsigned>(si.scany());
03652         for (int x = si.startx(); x<=si.endx(); ++x)
03653         {
03654           unsigned i = static_cast<unsigned>(x);
03655           vil_image_view<unsigned short> v = img->get_view(i, 1,j, 1);
03656           float fv = static_cast<float>(v(0,0));
03657           if (fv<min) min = fv;
03658           if (fv>max) max = fv;
03659           pixels.push_back(fv);
03660         }
03661       }
03662       return pixels;
03663     }
03664     else
03665     {
03666       for (si.reset(); si.next();)
03667       {
03668         unsigned j = static_cast<unsigned>(si.scany());
03669         for (int x = si.startx(); x<=si.endx(); ++x)
03670         {
03671           unsigned i = static_cast<unsigned>(x);
03672           vil_image_view<unsigned short> v = img->get_view(i, 1,j, 1);
03673           float fv = 0;
03674           for (unsigned p = 0; p<np; ++p)
03675             fv += v(0,0,p);
03676           fv/=3;
03677           if (fv<min) min = fv;
03678           if (fv>max) max = fv;
03679           pixels.push_back(fv);
03680         }
03681       }
03682       return pixels;
03683     }
03684   }
03685   vcl_cerr << "In brip_vil_float_ops::scan_region() - unknown format\n";
03686   return pixels;
03687 }
03688 
03689 // It has been observed that color order is somewhat invariant to illumination
03690 // Compute an index based on color order. The tolerance determines if two
03691 // color bands are too close to determine order, i.e. they should be considered
03692 // equal instead
03693 // the two relations being considered  are <, > and =, so the relationship
03694 // graph looks like:
03695 // \verbatim
03696 //         G
03697 //       /   \.
03698 //   > < =   > < =
03699 //    /         \.
03700 //  R  - > < = -  B
03701 // \endverbatim
03702 // Thus, there are three graph edges with each of three possible labels or
03703 // 9 possible order codes. An easy coding scheme is to use the top 6 bits of
03704 // the byte output pixel. The relationship is encoded as states of bit pairs
03705 // \verbatim
03706 // Color relations  R*G  R*B  G*B    * indicates > < = (1,2,3)
03707 // Bit indices      7,6  5,4  3,2
03708 // \endverbatim
03709 //
03710 vil_image_view<vxl_byte> brip_vil_float_ops::
03711 color_order(vil_image_view<float> const& color_image, float eq_tol)
03712 {
03713   unsigned ni = color_image.ni(), nj = color_image.nj();
03714   unsigned np = color_image.nplanes();
03715   vil_image_view<vxl_byte> temp;
03716   if (np!=3)
03717     return temp; // improper call
03718   temp.set_size(ni, nj);
03719   vxl_byte rg_eq=192, rg_lt=128, rg_gt=64;
03720   vxl_byte rb_eq=48, rb_lt=32, rb_gt=16;
03721   vxl_byte gb_eq=12, gb_lt=8, gb_gt=4;
03722   for (unsigned j = 0; j<nj; ++j)
03723     for (unsigned i = 0; i<ni; ++i) {
03724       float r=color_image(i,j,0), g=color_image(i,j,1), b=color_image(i,j,2);
03725       vxl_byte rg, rb, gb;
03726       // rg code
03727       if (vcl_fabs(r-g)<eq_tol)
03728         rg = rg_eq;
03729       else if (r<g)
03730         rg = rg_lt;
03731       else rg = rg_gt;
03732       // rb code
03733       if (vcl_fabs(r-b)<eq_tol)
03734         rb = rb_eq;
03735       else if (r<b)
03736         rb = rb_lt;
03737       else rb = rb_gt;
03738       // gb code
03739       if (vcl_fabs(g-b)<eq_tol)
03740         gb = gb_eq;
03741       else if (g<b)
03742         gb = gb_lt;
03743       else gb = gb_gt;
03744       vxl_byte v = static_cast<vxl_byte>(rg|rb|gb); // bitwise or
03745       temp(i,j) = v;
03746     }
03747   return temp;
03748 }
03749 
03750 #if 0 // only used in commented-out part of brip_vil_float_ops::extrema()
03751 static double zs(double x)
03752 { if (x < 0.0001 && x > -0.0001) return 0; else return x; }
03753 #endif // 0
03754 
03755 static double brip_vil_rot_gauss(double x, double y,
03756                                  double sigma_x, double sigma_y, double theta)
03757 {
03758   double theta_rad = theta*vnl_math::pi_over_180;
03759   double s = vcl_sin(theta_rad), c = vcl_cos(theta_rad);
03760   double sxr = 1.0/sigma_x, syr = 1.0/sigma_y;
03761   double ax = (c*x + s*y)*sxr, ay = (-s*x + c*y)*syr;
03762   double ret = vcl_exp(-0.5*(ax*ax + ay*ay));
03763   return ret*sxr*syr/vnl_math::twopi;
03764 }
03765 
03766 // revert angle to the range [-90, 90]
03767 float brip_vil_float_ops::extrema_revert_angle(float angle)
03768 {
03769   if (angle>=0.0f) {
03770     if (angle>90.0f&&angle<=270.0f)
03771       angle -= 180.0f;
03772     else if (angle>270.0f&&angle<=360.0f)
03773       angle -= 360.0f;
03774   }
03775   if (angle<0.0f) {
03776     if (angle<-90.0f&&angle>=-270.0f)
03777       angle += 180.0f;
03778     else if (angle<-270.0f&&angle>=-360.0f)
03779       angle += 360.0f;
03780   }
03781   return angle;
03782 }
03783 
03784 // if scale_invariant = true then the response mag is independent of lambda1
03785 void brip_vil_float_ops::extrema_kernel_mask(float lambda0, float lambda1,
03786                                              float theta,
03787                                              vbl_array_2d<float>& kernel,
03788                                              vbl_array_2d<bool>& mask, float cutoff_percentage, bool scale_invariant)
03789 {
03790   theta = extrema_revert_angle(theta); // in range [-90, 90]
03791   // convert theta to radians
03792   double theta_rad = theta*vnl_math::pi_over_180;
03793   double s = vcl_sin(theta_rad), c = vcl_cos(theta_rad);
03794   double s0 = lambda0, s1 = lambda1;
03795   double s1sq = 1.0/(s1*s1);
03796   // determine the size of the array
03797   double max_v = brip_vil_rot_gauss(0, 0, s0, s1, 0);
03798   double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
03799   unsigned ru = 0;
03800   while (brip_vil_rot_gauss(ru, 0, s0, s1, 0) >= cutoff)
03801     ++ru;
03802   unsigned rv = 0;
03803   while (brip_vil_rot_gauss(0, rv, s0, s1, 0) >= cutoff)
03804     ++rv;
03805 
03806   // rotate to get bounds
03807   int ri = static_cast<int>(vcl_fabs(ru*c+rv*s) +0.5);
03808   int rj = static_cast<int>(vcl_fabs(ru*s+rv*c) +0.5);
03809   if (s<0) {
03810     ri = static_cast<int>(vcl_fabs(ru*c-rv*s) +0.5);
03811     rj = static_cast<int>(vcl_fabs(ru*s-rv*c) +0.5);
03812   }
03813 
03814   unsigned ncols = 2*ri +1, nrows = 2*rj+1;
03815   vbl_array_2d<double> coef(nrows, ncols);
03816   mask.resize(nrows,ncols);
03817   double residual = 0.0, total  = 0.0;
03818   for (int ry = -rj; ry<=rj; ++ry)
03819     for (int rx = -ri; rx<=ri; ++rx)
03820     {
03821       double g = brip_vil_rot_gauss(rx, ry, s0, s1, theta);
03822       mask[ry+rj][rx+ri] = g>=cutoff;
03823       double temp = (-s*rx + c*ry);
03824       temp = temp*temp*s1sq;
03825       double v = (temp -1)*g*(scale_invariant?1.0:s1sq);
03826       if (g<cutoff) v = 0.0;
03827       coef[ry+rj][rx+ri] = v;
03828       residual+=v;
03829       total += vcl_fabs(v);
03830     }
03831   double cor = 0.0;
03832   if (total != 0)
03833     cor = -residual/total;
03834   kernel.resize(nrows, ncols);
03835   // correct any residual offset in coefficients
03836   // distribute proportionally to coefficient magnitude
03837   for (unsigned j = 0; j<nrows; ++j)
03838     for (unsigned i = 0; i<ncols; ++i)
03839     {
03840       double v = vcl_fabs(coef[j][i]);
03841       coef[j][i]+=v*cor;
03842       kernel[j][i]=static_cast<float>(coef[j][i]);
03843     }
03844 }
03845 
03846 //:
03847 //  \p theta must be given in degrees.
03848 //  Scale invariant means that the response is independent of the
03849 //  \a sigma_y of the unrotated operator, i.e. the direction of the derivative
03850 void brip_vil_float_ops::gaussian_kernel_mask(float lambda, vbl_array_2d<float>& kernel, vbl_array_2d<bool>& mask, float cutoff_percentage)
03851 {
03852   double s1 = lambda;
03853   double s1sq = 1.0/(s1*s1);
03854   // determine the size of the array
03855   double max_v = brip_vil_rot_gauss(0, 0, s1, s1, 0);
03856   double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
03857   unsigned r = 0;
03858   while (brip_vil_rot_gauss(r, 0, s1, s1, 0) >= cutoff)
03859     ++r;
03860 
03861   int ri = static_cast<int>(vcl_fabs((float)r) +0.5);
03862 
03863   unsigned n = 2*ri +1;
03864   mask.resize(n,n);
03865   kernel.resize(n, n);
03866   for (int ry = -ri; ry<=ri; ++ry)
03867     for (int rx = -ri; rx<=ri; ++rx)
03868     {
03869       double g = brip_vil_rot_gauss(rx, ry, s1, s1, 0);
03870       mask[ry+ri][rx+ri] = g>=cutoff;
03871       double temp = ry*s1sq;
03872       double v = (temp -1)*g;
03873       if (g<cutoff) v = 0.0;
03874       kernel[ry+ri][rx+ri] = static_cast<float>(vcl_fabs(v));
03875     }
03876 }
03877 
03878 //: return a square mask and the coefficient matrix [kernel] for a symmetric gaussian distribution
03879 void brip_vil_float_ops::gaussian_kernel_square_mask(float lambda, vbl_array_2d<float>& kernel, vbl_array_2d<bool>& mask, float cutoff_percentage)
03880 {
03881   double s1 = lambda;
03882   double s1sq = 1.0/(s1*s1);
03883   // determine the size of the array
03884   double max_v = brip_vil_rot_gauss(0, 0, s1, s1, 0);
03885   double cutoff = max_v*cutoff_percentage; // if 0.01 then 1% tails removed
03886   unsigned r = 0;
03887   while (brip_vil_rot_gauss(r, 0, s1, s1, 0) >= cutoff)
03888     ++r;
03889 
03890   int ri = static_cast<int>(vcl_fabs((float)r) +0.5);
03891 
03892   unsigned n = 2*ri +1;
03893   vbl_array_2d<double> coef(n, n);
03894   mask.resize(n,n);
03895   kernel.resize(n, n);
03896   for (int ry = -ri; ry<=ri; ++ry)
03897     for (int rx = -ri; rx<=ri; ++rx)
03898     {
03899       double g = brip_vil_rot_gauss(rx, ry, s1, s1, 0);
03900       mask[ry+ri][rx+ri] = true;
03901       double temp = ry*s1sq;
03902       kernel[ry+ri][rx+ri] = static_cast<float>(vcl_fabs((temp -1)*g));
03903     }
03904 }
03905 
03906 
03907 // Compute the standard deviation of an operator response
03908 // given the image intensity standard deviation at each pixel
03909 vil_image_view<float> brip_vil_float_ops::
03910 std_dev_operator(vil_image_view<float> const& sd_image,
03911                  vbl_array_2d<float> const& kernel)
03912 {
03913   unsigned ni = sd_image.ni(), nj = sd_image.nj();
03914   unsigned nrows = kernel.rows(), ncols = kernel.cols();
03915   int rrad = (nrows-1)/2, crad = (ncols-1)/2;
03916   float sd_min, sd_max;
03917   vil_math_value_range(sd_image, sd_min, sd_max);
03918 
03919   vil_image_view<float> res(ni, nj);
03920   res.fill(sd_max);
03921 
03922   vbl_array_2d<float>ksq(nrows, ncols);
03923   for (unsigned r = 0; r<nrows; ++r)
03924     for (unsigned c = 0; c<ncols; ++c)
03925       ksq[r][c] = kernel[r][c]*kernel[r][c];
03926 
03927   for (int j = rrad; j<static_cast<int>(nj-rrad); ++j)
03928     for (int i = crad; i<static_cast<int>(ni-crad); ++i)
03929     {
03930       float sum = 0;
03931       for (int jj = -rrad; jj<=rrad; ++jj)
03932         for (int ii = -crad; ii<=crad; ++ii) {
03933           float sd_sq = sd_image(i+ii, j+jj);
03934           sd_sq *= sd_sq;
03935           sum += sd_sq*ksq[jj+rrad][ii+crad];
03936         }
03937       res(i,j) = vcl_sqrt(sum);
03938     }
03939   return res;
03940 }
03941 
03942 // Compute the standard deviation of an operator response
03943 // given the image intensity standard deviation at each pixel
03944 // uses a modified formula to compute std_dev
03945 vil_image_view<float> brip_vil_float_ops::
03946 std_dev_operator_method2(vil_image_view<float> const& sd_image,
03947                          vbl_array_2d<float> const& kernel)
03948 {
03949   unsigned ni = sd_image.ni(), nj = sd_image.nj();
03950   unsigned nrows = kernel.rows(), ncols = kernel.cols();
03951   int rrad = (nrows-1)/2, crad = (ncols-1)/2;
03952   float sd_min, sd_max;
03953   vil_math_value_range(sd_image, sd_min, sd_max);
03954 
03955   vil_image_view<float> res(ni, nj);
03956   res.fill(sd_max);
03957 
03958   for (int j = rrad; j<static_cast<int>(nj-rrad); ++j)
03959     for (int i = crad; i<static_cast<int>(ni-crad); ++i)
03960     {
03961       float sum = 0;
03962       for (int jj = -rrad; jj<=rrad; ++jj)
03963         for (int ii = -crad; ii<=crad; ++ii) {
03964           float sd = sd_image(i+ii, j+jj);
03965           sum += (float)(sd*vcl_abs(kernel[jj+rrad][ii+crad]));
03966         }
03967       res(i,j) = sum;
03968     }
03969   return res;
03970 }
03971 
03972 vil_image_view<float>
03973 brip_vil_float_ops::extrema(vil_image_view<float> const& input,
03974                             float lambda0, float lambda1, float theta,
03975                             bool bright, bool mag_only,
03976                             bool output_response_mask,
03977                             bool signed_response, bool scale_invariant,
03978                             bool non_max_suppress, float cutoff_per)
03979 {
03980   assert(!(mag_only&&signed_response));
03981   vbl_array_2d<float> fa;
03982   vbl_array_2d<bool> mask;
03983   brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta,
03984                                           fa, mask, cutoff_per,
03985                                           scale_invariant);
03986   unsigned nrows = fa.rows(), ncols = fa.cols();
03987   int rj = (nrows-1)/2, ri = (ncols-1)/2;
03988   vbl_array_2d<double> coef(nrows,ncols);
03989   for (unsigned r = 0; r<nrows; ++r)
03990     for (unsigned c = 0; c<ncols; ++c)
03991       coef[r][c]=fa[r][c];
03992 
03993 #if 0
03994   // (unused) double max_v = brip_vil_rot_gauss(0, 0, lambda0, lambda1, 0);
03995   // (unused) double cutoff=static_cast<float>(max_v*0.01); // 1% tails removed
03996   vcl_cout << "\ngauss ={";
03997   for (unsigned j = 0; j<nrows; ++j) {
03998     vcl_cout << '{';
03999     for (unsigned i = 0; i<ncols-1; ++i)
04000       vcl_cout << zs(coef[j][i]) << ',';
04001     if (j != nrows-1)
04002       vcl_cout << zs(coef[j][ncols-1]) << "},";
04003     else
04004       vcl_cout << zs(coef[j][ncols-1]) << '}';
04005   }
04006 
04007   vcl_cout << "};\n\nmask ={";
04008   for (unsigned j = 0; j<nrows; ++j) {
04009     vcl_cout << '{';
04010     for (unsigned i = 0; i<ncols-1; ++i)
04011       vcl_cout << mask[j][i] << ',';
04012     if (j != nrows-1)
04013       vcl_cout << mask[j][ncols-1] << "},";
04014     else
04015       vcl_cout << mask[j][ncols-1] << '}';
04016   }
04017   vcl_cout << "};" << vcl_flush;
04018 #endif
04019   // compute response image
04020   unsigned ni = input.ni(), nj = input.nj();
04021   vil_image_view<float> temp(ni, nj);
04022   vil_image_view<float> temp2(ni, nj);
04023   temp.fill(0.0f); temp2.fill(0.0f);
04024   for (unsigned j = rj; j<(nj-rj); j++)
04025     for (unsigned i = ri; i<(ni-ri); i++) {
04026       double sum = 0;
04027       for (int jj=-rj; jj<=rj; ++jj)
04028         for (int ii=-ri; ii<=ri; ++ii)
04029           if (mask[jj+rj][ii+ri])
04030             sum += coef[jj+rj][ii+ri]*input(i+ii, j+jj);
04031       temp2(i,j) = static_cast<float>(sum);
04032       if (mag_only) {
04033         temp(i,j) = static_cast<float>(vcl_fabs(sum));
04034       }
04035       else if (bright) { // coefficients are negative at center
04036         if (sum<0) temp(i,j) = static_cast<float>(-sum);
04037       }
04038       else {
04039         if (sum>0) temp(i,j) = static_cast<float>(sum);
04040       }
04041     }
04042   // non-max suppression
04043   vil_image_view<float> res(temp);
04044   if (non_max_suppress)
04045     for (unsigned j = rj; j<(nj-rj); j++)
04046       for (unsigned i = ri; i<(ni-ri); i++)
04047       {
04048         float cv = temp(i,j);
04049         for (int jj=-rj; jj<=rj; ++jj)
04050           for (int ii=-ri; ii<=ri; ++ii)
04051             if ((ii!=0||jj!=0) && mask[jj+rj][ii+ri] && temp(i+ii, j+jj)>cv)
04052               res(i,j)=0.0f;
04053       }
04054   if (!output_response_mask&&!signed_response)
04055     return res;
04056   unsigned np = 2;
04057   if (output_response_mask&&signed_response)
04058     np = 3;
04059   //response plane and mask plane and/or signed response
04060   vil_image_view<float> res_mask(ni, nj, np);
04061   res_mask.fill(0.0f);
04062   // always copy response to plane 0
04063   for (unsigned j = rj; j<(nj-rj); j++)
04064     for (unsigned i = ri; i<(ni-ri); i++)
04065     {
04066       float rv = res(i,j);
04067       res_mask(i,j,0)=rv;
04068     }
04069   // copy mask plane to plane 1 (or not)
04070   if (output_response_mask) {
04071     for (unsigned j = rj; j<(nj-rj); j++)
04072       for (unsigned i = ri; i<(ni-ri); i++)
04073       {
04074         float rv = res(i,j);
04075         if (!rv)
04076           continue;
04077         for (int jj=-rj; jj<=rj; ++jj)
04078           for (int ii=-ri; ii<=ri; ++ii)
04079             if (mask[jj+rj][ii+ri])
04080               if (rv>res_mask(i+ii,j+jj,1))
04081                 res_mask(i+ii,j+jj,1) = rv;
04082       }
04083   }
04084   // copy signed response to plane 1 if no mask plane, otherwise plane 2
04085   if (signed_response) {
04086     unsigned p = np-1;
04087     for (unsigned j = rj; j<(nj-rj); j++)
04088       for (unsigned i = ri; i<(ni-ri); i++)
04089         res_mask(i,j,p) = temp2(i,j);
04090   }
04091   return res_mask;
04092 }
04093 
04094 
04095 //: Find anisotropic intensity extrema at a range of orientations and return the maximal response at the best orientation. Theta interval is in degrees
04096 //  if lambda0 == lambda1 then reduces to the normal extrema operator
04097 vil_image_view<float> brip_vil_float_ops::
04098 extrema_rotational(vil_image_view<float> const& input, float lambda0,
04099                    float lambda1, float theta_interval, bool bright,
04100                    bool mag_only, bool signed_response,
04101                    bool scale_invariant, bool non_max_suppress,
04102                    float cutoff_per)
04103 {
04104   //invalid response
04105   assert(!(mag_only&&signed_response));
04106 
04107   unsigned ni = input.ni();
04108   unsigned nj = input.nj();
04109 
04110   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
04111   output.fill(0.0f);
04112 
04113   if (lambda0 == lambda1) {  // if isotropic reduces to normal extrema calculation (perfect orientation symmetry)
04114     vil_image_view<float> res_mask_current =
04115       extrema(input, lambda0, lambda1, 0.0f, bright, mag_only,
04116               true, signed_response, scale_invariant,
04117               non_max_suppress, cutoff_per);
04118     unsigned np = res_mask_current.nplanes();
04119     for (unsigned j = 0; j<nj; j++)
04120       for (unsigned i = 0; i<ni; i++) {
04121         if (!signed_response)
04122           output(i,j,0) = res_mask_current(i,j,0);  // return the non-max suppressed for angle 0
04123         else
04124           output(i,j,0) = res_mask_current(i,j,np-1);
04125         output(i,j,1) = 0.0f;
04126         output(i,j,2) = res_mask_current(i,j,1);
04127       }
04128     return output;
04129   }
04130 
04131   // The kernel generator does not treat the x and y axis symmetrically, the method works correctly only when lambda0 > lambda1
04132   // 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
04133   if (lambda0 < lambda1) {
04134     vcl_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";
04135     throw 0;
04136   }
04137   if (theta_interval < vcl_numeric_limits<float>::epsilon()) {
04138     vcl_cout << "In brip_vil_float_ops::extrema_rotational() - ERROR! theta_interval needs to be larger than 0!\n";
04139     throw 0;
04140   }
04141 
04142   vil_image_view<float> res_img(ni, nj);
04143   vil_image_view<int> res_angle(ni, nj);
04144   res_img.fill(0.0f); res_angle.fill(0);
04145 
04146   vcl_vector<float> angles;
04147   angles.push_back(-1.0f);
04148   for (float theta = 0.0f; theta < 180.0f; theta += theta_interval) { angles.push_back(theta); }
04149 
04150   vcl_vector<vbl_array_2d<bool> > mask_vect(angles.size(), vbl_array_2d<bool>());
04151   int max_rji = 0;
04152   // elliptical operator has 180 degree rotational symmetry, so only the angles in the range [0,180] matter
04153   float theta = 0.0f; unsigned theta_i = 1;
04154   for ( ; theta < 180; theta += theta_interval, theta_i++)
04155   {
04156     // compute the response
04157     vbl_array_2d<float> fa; vbl_array_2d<bool> mask;
04158     brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta,
04159                                             fa, mask, cutoff_per, scale_invariant);
04160     mask_vect[theta_i] = mask;
04161     unsigned nrows = fa.rows(), ncols = fa.cols();
04162     int rj = (nrows-1)/2, ri = (ncols-1)/2;
04163     if (rj > max_rji) max_rji = rj;
04164     if (ri > max_rji) max_rji = ri;
04165     for (unsigned j = rj; j<(nj-rj); j++) {
04166       for (unsigned i = ri; i<(ni-ri); i++) {
04167         float res = 0.0f;
04168         double sum = 0;
04169         for (int jj=-rj; jj<=rj; ++jj)
04170           for (int ii=-ri; ii<=ri; ++ii)
04171             if (mask[jj+rj][ii+ri]) {
04172               sum += double(fa[jj+rj][ii+ri])*input(i+ii, j+jj);
04173             }
04174         if (mag_only) {
04175           res = static_cast<float>(vcl_fabs(sum));
04176         }
04177         else if (signed_response) {
04178           res = static_cast<float>(sum);
04179         }
04180         else if (bright) { // coefficients are negative at center
04181           if (sum<0) res = static_cast<float>(-sum);
04182         }
04183         else {
04184           if (sum>0) res = static_cast<float>(sum);
04185         }
04186         // pick largest angle response magnitude
04187         if (vcl_fabs(res_img(i,j)) < vcl_fabs(res)) {
04188           res_img(i,j) = res;
04189           res_angle(i,j) = theta_i;
04190         }
04191       }
04192     }
04193     vcl_cout << '.';
04194   }
04195   vcl_cout << '\n';
04196   if (!non_max_suppress) return res_img;
04197   // now we have pixel-wise best angle, run the non-max suppression around each non-zero pixel using the angles mask
04198   vil_image_view<float> res(res_img);
04199 
04200   for (unsigned j = max_rji; j<(nj-max_rji); j++)
04201     for (unsigned i = max_rji; i<(ni-max_rji); i++)
04202     {
04203       float cv = res_img(i,j);
04204       if (!cv)
04205         continue;
04206       int theta_i = res_angle(i,j);
04207       vbl_array_2d<bool> mask = mask_vect[theta_i];
04208       unsigned nrows = mask.rows(), ncols = mask.cols();
04209       int rj = (nrows-1)/2, ri = (ncols-1)/2;
04210 
04211       bool max = true;
04212       for (int jj=-rj; jj<=rj; ++jj) {
04213         for (int ii=-ri; ii<=ri; ++ii) {
04214           if ((ii!=0||jj!=0) && mask[jj+rj][ii+ri] && res_img(i+ii, j+jj)>cv) {
04215             max = false;
04216             break;
04217           }
04218         }
04219         if (!max)
04220           break;
04221       }
04222       if (!max) {
04223         res(i,j) = 0.0f;
04224         res_angle(i, j) = 0; // the zeroth angle is -1.0f, so invalid
04225       }
04226     }
04227 
04228   vil_image_view<float> res_mask(ni, nj);
04229   res_mask.fill(0.0f);
04230   // now all the non-zero elements in res are our responses, create the mask image using the angle info
04231   for (unsigned j = max_rji; j<(nj-max_rji); j++)
04232     for (unsigned i = max_rji; i<(ni-max_rji); i++)
04233     {
04234       float rv = res(i,j);
04235       if (!rv)
04236         continue;
04237       int theta_i = res_angle(i,j);
04238       // get the mask for this angle
04239       vbl_array_2d<bool> mask = mask_vect[theta_i];
04240       unsigned nrows = mask.rows(), ncols = mask.cols();
04241       int rj = (nrows-1)/2, ri = (ncols-1)/2;
04242 
04243       for (int jj=-rj; jj<=rj; ++jj)
04244         for (int ii=-ri; ii<=ri; ++ii)
04245           if (mask[jj+rj][ii+ri])
04246             if (rv>res_mask(i+ii,j+jj))
04247               res_mask(i+ii,j+jj) = rv;
04248     }
04249 
04250   // now prepare the output accordingly
04251   for (unsigned j = max_rji; j<(nj-max_rji); j++)
04252     for (unsigned i = max_rji; i<(ni-max_rji); i++)
04253     {
04254       output(i,j,0) = res(i,j);
04255       output(i,j,1) = angles[res_angle(i,j)];
04256       output(i,j,2) = res_mask(i,j);
04257     }
04258   return output;
04259 }
04260 
04261 
04262 //: theta and phi are in radians
04263 float brip_vil_float_ops::elu(float phi, float lamda0,
04264                               float lambda1, float theta)
04265 {
04266   double sp = vcl_sin(phi), cp = vcl_cos(phi);
04267   double st = vcl_sin(theta), ct = vcl_cos(theta);
04268   double temp = 3.0*lamda0*cp*ct-3.0*lambda1*sp*st;
04269   return static_cast<float>(temp);
04270 }
04271 
04272 //: theta and phi are in radians
04273 float brip_vil_float_ops::elv(float phi, float lamda0,
04274                               float lambda1, float theta)
04275 {
04276   double sp = vcl_sin(phi), cp = vcl_cos(phi);
04277   double st = vcl_sin(theta), ct = vcl_cos(theta);
04278   double temp = 3.0*lamda0*st*cp + 3.0*lambda1*sp*ct;
04279   return static_cast<float>(temp);
04280 }
04281 
04282 //: theta and phi are in radians
04283 void brip_vil_float_ops::
04284 max_inscribed_rect(float lambda0, float lambda1, float theta,
04285                    float& u_rect, float& v_rect)
04286 {
04287   float sign = -1.0f;
04288   if (theta<0) sign = 1.0f;
04289   float th_rad = static_cast<float>(theta*vnl_math::pi_over_180);
04290   float maxa = 0.0f, max_phi = 0.0f;
04291   for (float phi = -float(vnl_math::pi); phi<=float(vnl_math::pi); phi+=0.01f)
04292   {
04293     float a = (1.0f+brip_vil_float_ops::elu(phi, lambda0, lambda1, th_rad));
04294     a *= (1.0f+ sign*brip_vil_float_ops::elv(phi, lambda0, lambda1, th_rad));
04295     if (a>maxa)
04296     {
04297       maxa = a;
04298       max_phi = phi;
04299     }
04300   }
04301   u_rect = brip_vil_float_ops::elu(max_phi, lambda0, lambda1, th_rad);
04302   v_rect = brip_vil_float_ops::elv(max_phi, lambda0, lambda1, th_rad);
04303   v_rect = vcl_fabs(v_rect);
04304 }
04305 
04306 static void eu(float lambda0, float cutoff_per, vcl_vector<float>& kernel)
04307 {
04308   kernel.clear();
04309   double l_sqi = 1.0/(lambda0*lambda0);
04310   double norm = vnl_math::one_over_sqrt2pi;
04311   norm /= lambda0;
04312   int r = static_cast<int>(vcl_sqrt((-2.0*vcl_log(cutoff_per))/l_sqi)+0.5);
04313   for (int i = -r; i<=r; ++i)
04314   {
04315     double k = vcl_exp(-0.5*i*i*l_sqi);
04316     k*=norm;
04317     kernel.push_back(static_cast<float>(k));
04318   }
04319 }
04320 
04321 static void ev(float lambda1, float cutoff_per, bool scale_invariant,
04322                vcl_vector<float>& kernel)
04323 {
04324   kernel.clear();
04325   double l1_sqi = 1/(lambda1*lambda1);
04326   int r = static_cast<int>(vcl_sqrt((-2.0*vcl_log(cutoff_per))/l1_sqi)+0.5);
04327   double norm = scale_invariant ? 1/lambda1 : l1_sqi/lambda1;
04328   norm *= vnl_math::one_over_sqrt2pi;
04329   for (int i = -r; i<=r; ++i)
04330   {
04331     double s = i*i*l1_sqi;
04332     double k = vcl_exp(-0.5*s);
04333     k *= norm*(s-1.0);
04334     kernel.push_back(static_cast<float>(k));
04335   }
04336 }
04337 
04338 static vil_image_view<float> convolve_u(vil_image_view<float> const& image,
04339                                         vcl_vector<float> const& kernel,
04340                                         bool initialize = true)
04341 {
04342   int nk = kernel.size();
04343   if (nk<3)
04344     return image;
04345   int ni = image.ni(), nj = image.nj();
04346   int ru = (nk-1)/2;
04347   vil_image_view<float> Iu(ni, nj);
04348   if (initialize) Iu.fill(0.0f);
04349   for (int j = 0; j<nj; ++j)
04350     for (int i = ru; i<ni-ru; ++i) {
04351       float sum = 0.0f;
04352       for (int k = -ru; k<=ru; ++k)
04353       {
04354         float v = image(i+k, j);
04355         sum += kernel[k+ru]*v;
04356       }
04357       Iu(i,j) = sum;
04358     }
04359   return Iu;
04360 }
04361 
04362 static vil_image_view<float> convolve_v(vil_image_view<float> const& image,
04363                                         vcl_vector<float> const& kernel,
04364                                         bool initialize = true)
04365 {
04366   int nk = kernel.size();
04367   if (nk<3)
04368     return image;
04369   int ni = image.ni(), nj = image.nj();
04370   int rv = (nk-1)/2;
04371   vil_image_view<float> Iv(ni, nj);
04372   if (initialize) Iv.fill(0.0f);
04373   for (int j = rv; j<nj-rv; ++j)
04374     for (int i = 0; i<ni; ++i) {
04375       float sum = 0.0f;
04376       for (int k = -rv; k<=rv; ++k)
04377       {
04378         float v = image(i, j+k);
04379         sum += kernel[k+rv]*v;
04380       }
04381       Iv(i,j) = sum;
04382     }
04383   return Iv;
04384 }
04385 
04386 // this function tracks the image origin during rotation
04387 // the calculations are the same as used in the homography method.
04388 static void rotation_offset(int ni, int nj, float theta_deg,
04389                             int i, int j, int& ti, int& tj)
04390 {
04391   ti= 0; tj = 0;
04392   double deg_to_rad = vnl_math::pi_over_180;
04393   double rang = deg_to_rad*theta_deg;
04394   vgl_h_matrix_2d<double> H;
04395   H.set_identity().set_rotation(rang);
04396   vsol_box_2d_sptr input_roi;
04397   vsol_polygon_2d_sptr input_poly, output_poly;
04398   input_roi = new vsol_box_2d();
04399   input_roi->add_point(0, 0);
04400   input_roi->add_point(ni, nj);
04401   input_poly = bsol_algs::poly_from_box(input_roi);
04402   if (!bsol_algs::homography(input_poly, H, output_poly))
04403     return;
04404   vsol_box_2d_sptr temp = output_poly->get_bounding_box();
04405   vgl_h_matrix_2d<double> t; t.set_identity().set_translation(-temp->get_min_x(),-temp->get_min_y());
04406   vgl_homg_point_2d<double> torg = (t*H)*vgl_homg_point_2d<double>(double(i),double(j));
04407   ti = static_cast<int>(torg.x());
04408   tj = static_cast<int>(torg.y());
04409 }
04410 
04411 vil_image_view<float> brip_vil_float_ops::
04412 fast_extrema(vil_image_view<float> const& input, float lambda0, float lambda1,
04413              float theta, bool bright, bool mag_only,bool output_response_mask,
04414              bool signed_response, bool scale_invariant, bool non_max_suppress,
04415              float cutoff)
04416 {
04417   // invalid choice
04418   assert(!(mag_only&&signed_response));
04419   // the kernels for u and v
04420   vcl_vector<float> euk, evk;
04421   eu(lambda0, cutoff, euk);
04422   ev(lambda1, cutoff, scale_invariant, evk);
04423   vil_image_view<float> pre_res;
04424   int tuf=0, tvf=0, tui=0, tvi=0;
04425   vil_image_view<float> rot = input;
04426   bool ang90 = theta==90.0f||theta == 270.0f||
04427     theta ==-90.0f||theta == -270.0f;
04428   // rotate the image to the specified angle (minus due to flipped v axis)
04429   if (theta!=0&&!ang90)
04430     rot = brip_vil_float_ops::rotate(input, -theta);
04431   if (!ang90) {
04432     // convolve kernel across rows
04433     vil_image_view<float> rotu = convolve_u(rot, euk);
04434     // then convolve along columns
04435     vil_image_view<float> rotuv = convolve_v(rotu, evk);
04436     // rotate back
04437     if (theta!=0)
04438     {
04439       pre_res = brip_vil_float_ops::rotate(rotuv, +theta);
04440       // there is a border around the image to contain the rotated version
04441       // which should be removed.
04442       rotation_offset(input.ni(), input.nj(), -theta, 0, 0, tuf, tvf);
04443       rotation_offset(rotuv.ni(), rotuv.nj(), +theta, tuf, tvf, tui, tvi);
04444     }
04445     else pre_res = rotuv;
04446   }
04447   else {
04448     // convolve along columns
04449     vil_image_view<float> rotv = convolve_v(rot, euk);
04450     // convolve across rows
04451     vil_image_view<float> rotvu = convolve_u(rotv, evk);
04452     pre_res = rotvu;
04453   }
04454   int ni = input.ni(), nj = input.nj();
04455   vil_image_view<float> res(ni, nj);
04456   res.fill(0.0f);
04457   vil_image_view<float> temp2(ni, nj);
04458   temp2.fill(0.0f);
04459   for (int j = 0; j<nj; ++j)
04460     for (int i = 0; i<ni; ++i)
04461     {
04462       float v = pre_res(i + tui, j+ tvi);
04463       if (mag_only) {
04464         res(i,j) = vcl_fabs(v);
04465       }
04466       else if (bright) { // coefficients are negative at center
04467         if (v<0) res(i,j) = -v;
04468       }
04469       else {
04470         if (v>0) res(i,j) = v;
04471       }
04472       temp2(i,j) = v;// signed output
04473     }
04474   if (!non_max_suppress && !output_response_mask) {
04475     return signed_response ? temp2 : res;
04476   }
04477   // TO DO! handle non_max_suppress==false but response_mask==true
04478   // December 12, 2011 --- JLM
04479   //
04480   // the second derivative convolution is complete
04481   // non-maximal suppression
04482   // determine the inscribed rect in the response mask with max area
04483   float u_rect, v_rect;
04484   brip_vil_float_ops::max_inscribed_rect(lambda0,lambda1,theta,u_rect,v_rect);
04485   int ni_rect = static_cast<int>(u_rect);
04486   int nj_rect = static_cast<int>(v_rect);
04487   vil_image_view<float> res_sup(ni, nj);
04488   res_sup.fill(0.0f);
04489   for ( int i= ni_rect; i< ni-ni_rect; i+= ni_rect+1 )
04490     for ( int j = nj_rect; j < nj-nj_rect; j += nj_rect+1 ) {
04491       int ci= i, cj = j;
04492       bool find_response = false;
04493       for ( int di= 0; di<= ni_rect; di++ )
04494         for ( int dj = 0; dj <= nj_rect; dj++ )
04495         {
04496           float v = res(ci,cj);
04497           float dv = res(i+di,j+dj);
04498           if (v==0.0f&&dv==0.0f) continue;
04499           find_response = true;
04500           if ( v < dv )
04501           {
04502             ci= i+di;
04503             cj = j+dj;
04504           }
04505         }
04506       if (!find_response)
04507         continue;
04508       res_sup(ci,cj) = res(ci,cj);
04509     }
04510 
04511   // determine the outside bounding box
04512   vbl_array_2d<float> kern;
04513   vbl_array_2d<bool> mask;
04514   brip_vil_float_ops::extrema_kernel_mask(lambda0, lambda1, theta, kern, mask);
04515   int ri = (kern.cols()-1)/2, rj = (kern.rows()-1)/2;
04516   // go through the array again and search about each retained
04517   // local maximum using the outside bounding box
04518   // there will be at most one non-zero response to test in the inscribed rect
04519   for (int j = rj; j<(nj-rj); j++)
04520     for (int i = ri; i<(ni-ri); i++)
04521     {
04522       float v = res_sup(i,j);
04523       if (v==0.0f) continue;
04524       for (int jj=-rj; jj<=rj; ++jj)
04525         for (int ii=-ri; ii<=ri; ++ii)
04526           if (mask[rj+jj][ri+ii] && res_sup(i+ii, j+jj)<v)
04527             res_sup(i+ii,j+jj)=0.0f;
04528     }
04529 
04530   if (!output_response_mask&&!signed_response)
04531     return res_sup;
04532   unsigned np = 2;
04533   if (output_response_mask&&signed_response)
04534     np = 3;
04535   //response plane and mask plane and/or signed response
04536 
04537   vil_image_view<float> res_mask(ni, nj, np);
04538   res_mask.fill(0.0f);
04539   // always copy response to plane 0
04540   for (int j = rj; j<(nj-rj); j++)
04541     for (int i = ri; i<(ni-ri); i++)
04542     {
04543       float rv = res_sup(i,j);
04544       res_mask(i,j,0)=rv;
04545     }
04546   // copy mask plane to plane 1 (or not)
04547   if (output_response_mask) {
04548     for (int j = rj; j<(nj-rj); j++)
04549       for (int i = ri; i<(ni-ri); i++)
04550       {
04551         float rv = res_sup(i,j);
04552         if (!rv)
04553           continue;
04554         for (int jj=-rj; jj<=rj; ++jj)
04555           for (int ii=-ri; ii<=ri; ++ii)
04556             if (mask[jj+rj][ii+ri])
04557               if (rv>res_mask(i+ii,j+jj,1))
04558                 res_mask(i+ii,j+jj,1) = rv;
04559       }
04560   }
04561   // copy signed response to plane 1 if no mask plane, otherwise plane 2
04562   if (signed_response) {
04563     unsigned p = np-1;
04564     for (int j = rj; j<(nj-rj); j++)
04565       for (int i = ri; i<(ni-ri); i++)
04566         res_mask(i,j,p) = temp2(i,j);
04567   }
04568   return res_mask;
04569 }
04570 
04571 vil_image_view<float> brip_vil_float_ops::
04572 fast_extrema_rotational(vil_image_view<float> const& input,
04573                         float lambda0, float lambda1,
04574                         float theta_interval,
04575                         bool bright, bool mag_only,
04576                         bool signed_response, bool scale_invariant,
04577                         bool non_max_suppress,
04578                         float cutoff) {
04579   unsigned ni = input.ni(), nj = input.nj();
04580   vil_image_view<float> resp(ni, nj);
04581   resp.fill(0.0f);
04582 
04583     // elliptical operator has 180 degree rotational symmetry, so only the angles in the range [0,180] matter
04584   float theta = 0; unsigned theta_i = 1;
04585   for (; theta < 180; theta += theta_interval, theta_i++)
04586   {
04587    vil_image_view<float> temp =
04588      fast_extrema(input, lambda0, lambda1, theta, bright,
04589                   mag_only, false, signed_response,
04590                   scale_invariant, false, cutoff);
04591     for (unsigned j = 0; j<nj; ++j)
04592       for (unsigned i = 0; i<ni; ++i) {
04593         float v = temp(i,j);
04594         if (vcl_fabs(v)>vcl_fabs(resp(i,j)))
04595           resp(i,j) = v;
04596       }
04597   }
04598   return resp;
04599 }