[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [47622] branches/soc-2011-tomato/extern/ libmv: Make planar tracking much faster.
Keir Mierle
mierle at gmail.com
Fri Jun 8 19:42:25 CEST 2012
Revision: 47622
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=47622
Author: keir
Date: 2012-06-08 17:42:17 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
Make planar tracking much faster.
- This makes planar tracking around 2-3x or more faster than
before, by rearranging how the sampling is done.
Previously, the source patch was sampled repeatedly on
every optimizer iteration; this was done for
implementation speed, but was wasteful in computation.
- This also contains some additions to Ceres to help
deailing with mixed numeric / automatic differentation. In
particular, there is now a "Chain::Rule" operator that
facilitates calling a function that takes Jet arguments,
yet does numeric derivatives internally. This is used to
mix the numeric differentation of the images with the warp
parameters, passed as jets by Ceres to the warp functor.
There is also a new "JetOps" object for doing operations
on types which may or may not be jets, such as scaling
the derivative part only, or extracting the scalar part
of a jet.
The Ceres patches are aimed at upstream.
- A new function for sampling a patch is now part of the
track_region.h API; this will get used to make the preview
widget properly show what is getting tracked. Currently
the preview widget does not handle perspective tracks.
Known issues:
This patch introduces a bug such that the "Minimum
Correlation" flag does not work; if it is enabled, tracking
aborts immediately. The workaround for now is to disable the
correlation checking, and examine your tracks carefully. A
fix will get added shortly.
Modified Paths:
--------------
branches/soc-2011-tomato/extern/libmv/libmv/image/sample.h
branches/soc-2011-tomato/extern/libmv/libmv/tracking/track_region.cc
branches/soc-2011-tomato/extern/libmv/libmv/tracking/track_region.h
branches/soc-2011-tomato/extern/libmv/third_party/ceres/include/ceres/jet.h
Modified: branches/soc-2011-tomato/extern/libmv/libmv/image/sample.h
===================================================================
--- branches/soc-2011-tomato/extern/libmv/libmv/image/sample.h 2012-06-08 17:16:32 UTC (rev 47621)
+++ branches/soc-2011-tomato/extern/libmv/libmv/image/sample.h 2012-06-08 17:42:17 UTC (rev 47622)
@@ -37,15 +37,15 @@
static inline void LinearInitAxis(float fx, int width,
int *x1, int *x2,
float *dx1, float *dx2) {
- const int ix = int(fx);
+ const int ix = static_cast<int>(fx);
if (ix < 0) {
*x1 = 0;
*x2 = 0;
*dx1 = 1;
*dx2 = 0;
- } else if (ix > width-2) {
- *x1 = width-1;
- *x2 = width-1;
+ } else if (ix > width - 2) {
+ *x1 = width - 1;
+ *x2 = width - 1;
*dx1 = 1;
*dx2 = 0;
} else {
@@ -74,6 +74,27 @@
dy2 * ( dx1 * im21 + dx2 * im22 ));
}
+/// Linear interpolation, of all channels. The sample is assumed to have the
+/// same size as the number of channels in image.
+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;
+
+ LinearInitAxis(y, image.Height(), &y1, &y2, &dy1, &dy2);
+ LinearInitAxis(x, image.Width(), &x1, &x2, &dx1, &dx2);
+
+ for (int i = 0; i < image.Depth(); ++i) {
+ const T im11 = image(y1, x1, i);
+ const T im12 = image(y1, x2, i);
+ 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 ));
+ }
+}
+
// Downsample all channels by 2. If the image has odd width or height, the last
// row or column is ignored.
// FIXME(MatthiasF): this implementation shouldn't be in an interface file
Modified: branches/soc-2011-tomato/extern/libmv/libmv/tracking/track_region.cc
===================================================================
--- branches/soc-2011-tomato/extern/libmv/libmv/tracking/track_region.cc 2012-06-08 17:16:32 UTC (rev 47621)
+++ branches/soc-2011-tomato/extern/libmv/libmv/tracking/track_region.cc 2012-06-08 17:42:17 UTC (rev 47622)
@@ -43,6 +43,10 @@
namespace libmv {
+using ceres::Jet;
+using ceres::JetOps;
+using ceres::Chain;
+
TrackRegionOptions::TrackRegionOptions()
: mode(TRANSLATION),
minimum_correlation(0),
@@ -77,71 +81,31 @@
return true;
}
-// Because C++03 doesn't support partial template specializations for
-// functions, but at the same time member function specializations are not
-// supported, the sample function must be inside a template-specialized class
-// with a non-templated static member.
-
// 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.
-
-// Empty default template.
template<typename T>
-struct AutoDiff {
- // Sample only the image when the coordinates are scalars.
- static T Sample(const FloatImage &image_and_gradient,
- const T &x, const T &y) {
- return SampleLinear(image_and_gradient, y, x, 0);
- }
+static T SampleWithDerivative(const FloatImage &image_and_gradient,
+ const T &x,
+ const T &y) {
+ float scalar_x = JetOps<T>::GetScalar(x);
+ float scalar_y = JetOps<T>::GetScalar(y);
- static void SetScalarPart(double scalar, T *value) {
- *value = scalar;
+ // Note that sample[1] and sample[2] will be uninitialized in the scalar
+ // case, but that is not an issue because the Chain::Rule below will not read
+ // the uninitialized values.
+ float sample[3];
+ if (JetOps<T>::IsScalar()) {
+ // For the scalar case, only sample the image.
+ sample[0] = SampleLinear(image_and_gradient, scalar_y, scalar_x, 0);
+ } else {
+ // For the derivative case, sample the gradient as well.
+ SampleLinear(image_and_gradient, scalar_y, scalar_x, sample);
}
- static void ScaleDerivative(double scale_by, T *value) {
- // For double, there is no derivative to scale.
- }
-};
+ T xy[2] = { x, y };
+ return Chain<float, 2, T>::Rule(sample[0], sample + 1, xy);
+}
-// Sample the image and gradient when the coordinates are jets, applying the
-// jacobian appropriately to propagate the derivatives from the coordinates.
-template<typename T, int N>
-struct AutoDiff<ceres::Jet<T, N> > {
- static ceres::Jet<T, N> Sample(const FloatImage &image_and_gradient,
- const ceres::Jet<T, N> &x,
- const ceres::Jet<T, N> &y) {
- // Sample the image and its derivatives in x and y. One way to think of
- // this is that the image is a scalar function with a single vector
- // argument, xy, of dimension 2. Call this s(xy).
- const T s = SampleLinear(image_and_gradient, y.a, x.a, 0);
- const T dsdx = SampleLinear(image_and_gradient, y.a, x.a, 1);
- const T dsdy = SampleLinear(image_and_gradient, y.a, x.a, 2);
-
- // However, xy is itself a function of another variable ("z"); xy(z) =
- // [x(z), y(z)]^T. What this function needs to return is "s", but with the
- // derivative with respect to z attached to the jet. So combine the
- // derivative part of x and y's jets to form a Jacobian matrix between x, y
- // and z (i.e. dxy/dz).
- Eigen::Matrix<T, 2, N> dxydz;
- dxydz.row(0) = x.v.transpose();
- dxydz.row(1) = y.v.transpose();
-
- // Now apply the chain rule to obtain ds/dz. Combine the derivative with
- // the scalar part to obtain s with full derivative information.
- ceres::Jet<T, N> jet_s;
- jet_s.a = s;
- jet_s.v = Matrix<T, 1, 2>(dsdx, dsdy) * dxydz;
- return jet_s;
- }
-
- static void SetScalarPart(double scalar, ceres::Jet<T, N> *value) {
- value->a = scalar;
- }
- static void ScaleDerivative(double scale_by, ceres::Jet<T, N> *value) {
- value->v *= scale_by;
- }
-};
-
template<typename Warp>
class BoundaryCheckingCallback : public ceres::IterationCallback {
public:
@@ -188,36 +152,73 @@
canonical_to_image1_(canonical_to_image1),
num_samples_x_(num_samples_x),
num_samples_y_(num_samples_y),
- warp_(warp) {}
+ warp_(warp),
+ pattern_and_gradient_(num_samples_y_, num_samples_x_, 3),
+ pattern_positions_(num_samples_y_, num_samples_x_, 2),
+ pattern_mask_(num_samples_y_, num_samples_x_, 1) {
+ ComputeCanonicalPatchAndNormalizer();
+ }
+ void ComputeCanonicalPatchAndNormalizer() {
+ src_mean_ = 0.0;
+ double num_samples = 0.0;
+ for (int r = 0; r < num_samples_y_; ++r) {
+ for (int c = 0; c < num_samples_x_; ++c) {
+ // Compute the position; cache it.
+ Vec3 image_position = canonical_to_image1_ * Vec3(c, r, 1);
+ image_position /= image_position(2);
+ pattern_positions_(r, c, 0) = image_position(0);
+ pattern_positions_(r, c, 1) = image_position(1);
+
+ // Sample the pattern and gradients.
+ SampleLinear(image_and_gradient1_,
+ image_position(1), // Sample is r, c.
+ image_position(0),
+ &pattern_and_gradient_(r, c, 0));
+
+ // Sample sample the mask.
+ double mask_value = 1.0;
+ if (options_.image1_mask != NULL) {
+ SampleLinear(*options_.image1_mask,
+ image_position(1),
+ image_position(0),
+ &pattern_mask_(r, c, 0));
+ mask_value = pattern_mask_(r, c);
+ }
+ src_mean_ += pattern_and_gradient_(r, c, 0) * mask_value;
+ num_samples += mask_value;
+ }
+ }
+ src_mean_ /= num_samples;
+ }
+
template<typename T>
bool operator()(const T *warp_parameters, T *residuals) const {
+ if (options_.image1_mask != NULL) {
+ VLOG(2) << "Using a mask.";
+ }
for (int i = 0; i < Warp::NUM_PARAMETERS; ++i) {
VLOG(2) << "warp_parameters[" << i << "]: " << warp_parameters[i];
}
- T src_mean = T(1.0);
T dst_mean = T(1.0);
if (options_.use_normalized_intensities) {
- ComputeNormalizingCoefficients(warp_parameters,
- &src_mean,
- &dst_mean);
+ ComputeNormalizingCoefficient(warp_parameters,
+ &dst_mean);
}
int cursor = 0;
for (int r = 0; r < num_samples_y_; ++r) {
for (int c = 0; c < num_samples_x_; ++c) {
- // Compute the location of the source pixel (via homography).
- Vec3 image1_position = canonical_to_image1_ * Vec3(c, r, 1);
- image1_position /= image1_position(2);
-
+ // Use the pre-computed image1 position.
+ Vec2 image1_position(pattern_positions_(r, c, 0),
+ pattern_positions_(r, c, 1));
+
// Sample the mask early; if it's zero, this pixel has no effect. This
// allows early bailout from the expensive sampling that happens below.
double mask_value = 1.0;
if (options_.image1_mask != NULL) {
- mask_value = AutoDiff<double>::Sample(*options_.image1_mask,
- image1_position[0],
- image1_position[1]);
+ mask_value = pattern_mask_(r, c);
if (mask_value == 0.0) {
residuals[cursor++] = T(0.0);
continue;
@@ -233,54 +234,61 @@
&image2_position[1]);
// Sample the destination, propagating derivatives.
- T dst_sample = AutoDiff<T>::Sample(image_and_gradient2_,
- image2_position[0],
- image2_position[1]);
+ T dst_sample = SampleWithDerivative(image_and_gradient2_,
+ image2_position[0],
+ image2_position[1]);
// Sample the source. This is made complicated by ESM mode.
T src_sample;
- if (options_.use_esm) {
+ if (0 && 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
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list