Planar tracker polish.
authorKeir Mierle <mierle@gmail.com>
Sat, 9 Jun 2012 06:55:21 +0000 (06:55 +0000)
committerKeir Mierle <mierle@gmail.com>
Sat, 9 Jun 2012 06:55:21 +0000 (06:55 +0000)
- Fixes the correlation checking code that was broken in the
  previous commit. The bug was a transpose error.
- Fixes a memory leak of the warp functor, found by Sameer.
- Various cleanups done at Sameer's suggestion.

Thanks to Sameer Agarwal for a code review.

extern/libmv/libmv/image/sample.h
extern/libmv/libmv/tracking/track_region.cc
extern/libmv/libmv/tracking/track_region.h

index 7bcbe5b04013bef3e4bac6a90bb3f69faf98a6a4..04a5748ea4412559e0808064390a2cca124bc710 100644 (file)
@@ -34,25 +34,22 @@ inline T SampleNearest(const Array3D<T> &image,
   return image(i, j, v);
 }
 
-static inline void LinearInitAxis(float fx, int width,
-                                  int *x1, int *x2,
-                                  float *dx1, float *dx2) {
-  const int ix = static_cast<int>(fx);
+inline void LinearInitAxis(float x, int size,
+                           int *x1, int *x2,
+                           float *dx) {
+  const int ix = static_cast<int>(x);
   if (ix < 0) {
     *x1 = 0;
     *x2 = 0;
-    *dx1 = 1;
-    *dx2 = 0;
-  } else if (ix > width - 2) {
-    *x1 = width - 1;
-    *x2 = width - 1;
-    *dx1 = 1;
-    *dx2 = 0;
+    *dx = 1.0;
+  } else if (ix > size - 2) {
+    *x1 = size - 1;
+    *x2 = size - 1;
+    *dx = 1.0;
   } else {
     *x1 = ix;
-    *x2 = *x1 + 1;
-    *dx1 = *x2 - fx;
-    *dx2 = 1 - *dx1;
+    *x2 = ix + 1;
+    *dx = *x2 - x;
   }
 }
 
@@ -60,18 +57,18 @@ static inline void LinearInitAxis(float fx, int width,
 template<typename T>
 inline T SampleLinear(const Array3D<T> &image, float y, float x, int v = 0) {
   int x1, y1, x2, y2;
-  float dx1, dy1, dx2, dy2;
+  float dx, dy;
 
-  LinearInitAxis(y, image.Height(), &y1, &y2, &dy1, &dy2);
-  LinearInitAxis(x, image.Width(),  &x1, &x2, &dx1, &dx2);
+  LinearInitAxis(y, image.Height(), &y1, &y2, &dy);
+  LinearInitAxis(x, image.Width(),  &x1, &x2, &dx);
 
   const T im11 = image(y1, x1, v);
   const T im12 = image(y1, x2, v);
   const T im21 = image(y2, x1, v);
   const T im22 = image(y2, x2, v);
 
-  return T(dy1 * ( dx1 * im11 + dx2 * im12 ) +
-           dy2 * ( dx1 * im21 + dx2 * im22 ));
+  return T(     dy  * ( dx * im11 + (1.0 - dx) * im12 ) +
+           (1 - dy) * ( dx * im21 + (1.0 - dx) * im22 ));
 }
 
 /// Linear interpolation, of all channels. The sample is assumed to have the
@@ -79,10 +76,10 @@ inline T SampleLinear(const Array3D<T> &image, float y, float x, int v = 0) {
 template<typename T>
 inline void SampleLinear(const Array3D<T> &image, float y, float x, T *sample) {
   int x1, y1, x2, y2;
-  float dx1, dy1, dx2, dy2;
+  float dx, dy;
 
-  LinearInitAxis(y, image.Height(), &y1, &y2, &dy1, &dy2);
-  LinearInitAxis(x, image.Width(),  &x1, &x2, &dx1, &dx2);
+  LinearInitAxis(y, image.Height(), &y1, &y2, &dy);
+  LinearInitAxis(x, image.Width(),  &x1, &x2, &dx);
 
   for (int i = 0; i < image.Depth(); ++i) {
     const T im11 = image(y1, x1, i);
@@ -90,8 +87,8 @@ inline void SampleLinear(const Array3D<T> &image, float y, float x, T *sample) {
     const T im21 = image(y2, x1, i);
     const T im22 = image(y2, x2, i);
 
-    sample[i] = T(dy1 * ( dx1 * im11 + dx2 * im12 ) +
-                  dy2 * ( dx1 * im21 + dx2 * im22 ));
+    sample[i] = T(     dy  * ( dx * im11 + (1.0 - dx) * im12 ) +
+                  (1 - dy) * ( dx * im21 + (1.0 - dx) * im22 ));
   }
 }
 
index bcfb976fc68d2bf453f6dadf3f5bf984de8ce5b3..b9e4f883e78e21002dc4ac0fc4a3ff45cde20cf2 100644 (file)
@@ -81,9 +81,9 @@ bool AllInBounds(const FloatImage &image,
   return true;
 }
 
-// The "AutoDiff::Sample()" function allows sampling an image at an x, y
-// position such that if x and y are jets, then the derivative information is
-// correctly propagated.
+// Sample the image at position (x, y) but use the gradient, if present, to
+// propagate derivatives from x and y. This is needed to integrate the numeric
+// image gradients with Ceres's autodiff framework.
 template<typename T>
 static T SampleWithDerivative(const FloatImage &image_and_gradient,
                               const T &x,
@@ -172,7 +172,7 @@ class WarpCostFunctor {
 
         // Sample the pattern and gradients.
         SampleLinear(image_and_gradient1_,
-                     image_position(1),  // Sample is r, c.
+                     image_position(1),  // SampleLinear is r, c.
                      image_position(0),
                      &pattern_and_gradient_(r, c, 0));
 
@@ -180,7 +180,7 @@ class WarpCostFunctor {
         double mask_value = 1.0;
         if (options_.image1_mask != NULL) {
           SampleLinear(*options_.image1_mask,
-                       image_position(1),
+                       image_position(1),  // SampleLinear is r, c.
                        image_position(0),
                        &pattern_mask_(r, c, 0));
           mask_value = pattern_mask_(r, c);
@@ -216,6 +216,16 @@ class WarpCostFunctor {
 
         // Sample the mask early; if it's zero, this pixel has no effect. This
         // allows early bailout from the expensive sampling that happens below.
+        //
+        // Note that partial masks are not short circuited. To see why short
+        // circuiting produces bitwise-exact same results, consider that the
+        // residual for each pixel is 
+        //
+        //    residual = mask * (src - dst)  ,
+        //
+        // and for jets, multiplying by a scalar multiplies the derivative
+        // components by the scalar as well. Therefore, if the mask is exactly
+        // zero, then so too will the final residual and derivatives.
         double mask_value = 1.0;
         if (options_.image1_mask != NULL) {
           mask_value = pattern_mask_(r, c);
@@ -240,7 +250,7 @@ class WarpCostFunctor {
 
         // Sample the source. This is made complicated by ESM mode.
         T src_sample;
-        if (0 && options_.use_esm && !JetOps<T>::IsScalar()) {
+        if (options_.use_esm && !JetOps<T>::IsScalar()) {
           // In ESM mode, the derivative of the source is also taken into
           // account. This changes the linearization in a way that causes
           // better convergence. Copy the derivative of the warp parameters
@@ -270,8 +280,6 @@ class WarpCostFunctor {
           src_sample = T(pattern_and_gradient_(r, c));
         }
 
-        //LG << "src_sample: " << src_sample;
-
         // Normalize the samples by the mean values of each signal. The typical
         // light model assumes multiplicative intensity changes with changing
         // light, so this is a reasonable choice. Note that dst_mean has
@@ -281,8 +289,6 @@ class WarpCostFunctor {
           dst_sample /= dst_mean;
         }
 
-        //LG << "dst_sample: " << dst_sample;
-
         // The difference is the error.
         T error = src_sample - dst_sample;
 
@@ -389,8 +395,8 @@ class WarpCostFunctor {
 
         double x = pattern_and_gradient_(r, c);
         double y = SampleLinear(image_and_gradient2_,
-                                image2_position[0],
-                                image2_position[1]);
+                                image2_position[1],  // SampleLinear is r, c.
+                                image2_position[0]);
 
         // Weight the signals by the mask, if one is present.
         if (options_.image1_mask != NULL) {
@@ -1111,13 +1117,14 @@ void TemplatedTrackRegion(const FloatImage &image1,
                                          x1, y1, x2, y2);
     for (int i = 0; i < 4; ++i) {
       LG << "P" << i << ": (" << x1[i] << ", " << y1[i] << "); brute ("
-         << x2[i] << ", " << y2[i] << "); (dx, dy): (" << (x2[i] - x1[i]) << ", "
-         << (y2[i] - y1[i]) << ").";
+         << x2[i] << ", " << y2[i] << "); (dx, dy): (" << (x2[i] - x1[i])
+         << ", " << (y2[i] - y1[i]) << ").";
     }
   }
 
   // Prepare the initial warp parameters from the four correspondences.
-  // Note: This must happen after the brute initialization runs.
+  // Note: This must happen after the brute initialization runs, since the
+  // brute initialization mutates x2 and y2 in place.
   Warp warp(x1, y1, x2, y2);
 
   // Decide how many samples to use in the x and y dimensions.
@@ -1125,22 +1132,6 @@ void TemplatedTrackRegion(const FloatImage &image1,
   int num_samples_y;
   PickSampling(x1, y1, x2, y2, &num_samples_x, &num_samples_y);
 
-  ceres::Solver::Options solver_options;
-  solver_options.linear_solver_type = ceres::DENSE_QR;
-  solver_options.max_num_iterations = options.max_iterations;
-  solver_options.update_state_every_iteration = true;
-  solver_options.parameter_tolerance = 1e-16;
-  solver_options.function_tolerance = 1e-16;
-
-  // TODO(keir): Consider removing these options before committing.
-  solver_options.numeric_derivative_relative_step_size = 1e-3;
-  solver_options.check_gradients = false;
-  solver_options.gradient_check_relative_precision = 1e-10;
-  solver_options.minimizer_progress_to_stdout = false;
-
-  // Prevent the corners from going outside the destination image.
-  BoundaryCheckingCallback<Warp> callback(image2, warp, x1, y1);
-  solver_options.callbacks.push_back(&callback);
 
   // Compute the warp from rectangular coordinates.
   Mat3 canonical_homography = ComputeCanonicalHomography(x1, y1,
@@ -1158,9 +1149,7 @@ void TemplatedTrackRegion(const FloatImage &image1,
                                 warp);
 
   // Construct the problem with a single residual.
-  ceres::Problem::Options problem_options;
-  problem_options.cost_function_ownership = ceres::DO_NOT_TAKE_OWNERSHIP;
-  ceres::Problem problem(problem_options);
+  ceres::Problem problem;
   problem.AddResidualBlock(
       new ceres::AutoDiffCostFunction<
           WarpCostFunctor<Warp>,
@@ -1170,6 +1159,19 @@ void TemplatedTrackRegion(const FloatImage &image1,
       NULL,
       warp.parameters);
 
+  // Configure the solve.
+  ceres::Solver::Options solver_options;
+  solver_options.linear_solver_type = ceres::DENSE_QR;
+  solver_options.max_num_iterations = options.max_iterations;
+  solver_options.update_state_every_iteration = true;
+  solver_options.parameter_tolerance = 1e-16;
+  solver_options.function_tolerance = 1e-16;
+
+  // Prevent the corners from going outside the destination image.
+  BoundaryCheckingCallback<Warp> callback(image2, warp, x1, y1);
+  solver_options.callbacks.push_back(&callback);
+
+  // Run the solve.
   ceres::Solver::Summary summary;
   ceres::Solve(solver_options, &problem, &summary);
 
index ca21090e7c096d28d9b31cb7f84319796572c5ae..c23cde6000fad861779644b5e86480242a7f274c 100644 (file)
@@ -114,10 +114,10 @@ void TrackRegion(const FloatImage &image1,
                  TrackRegionResult *result);
 
 // Sample a "canonical" version of the passed planar patch, using bilinear
-// sampling. The passed corners must be within the image, possibly with a small
-// amount of slop, perhaps 2 pixels, around the edges (so e.g. a corner of the
-// patch cannot lie directly on the edge of the image). Four corners are always
-// required. All channels are interpolated.
+// sampling. The passed corners must be within the image, and have at least two
+// pixels of border around them. (so e.g. a corner of the patch cannot lie
+// directly on the edge of the image). Four corners are always required. All
+// channels are interpolated.
 bool SamplePlanarPatch(const FloatImage &image,
                        const double *xs, const double *ys,
                        int num_samples_x, int num_samples_y,