1 files changed, 31 insertions, 15 deletions

diff --git a/extern/libmv/libmv/multiview/euclidean_resection.cc b/extern/libmv/libmv/multiview/euclidean_resection.cc
index 92862515d7e..2605bf04622 100644
--- a/extern/libmv/libmv/multiview/euclidean_resection.cc
+++ b/extern/libmv/libmv/multiview/euclidean_resection.cc
@@ -37,13 +37,14 @@ typedef unsigned int uint;
 bool EuclideanResection(const Mat2X &x_camera, 
                         const Mat3X &X_world,
                         Mat3 *R, Vec3 *t,
-                        ResectionMethod method) {
+                        ResectionMethod method,
+                        double success_threshold) {
   switch (method) {
     case RESECTION_ANSAR_DANIILIDIS:
       EuclideanResectionAnsarDaniilidis(x_camera, X_world, R, t);
       break;
     case RESECTION_EPNP:
-      return EuclideanResectionEPnP(x_camera, X_world, R, t);      
+      return EuclideanResectionEPnP(x_camera, X_world, R, t, success_threshold);
       break;
     default:
       LOG(FATAL) << "Unknown resection method.";
@@ -351,9 +352,9 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
 }
 
 // Selects 4 virtual control points using mean and PCA.
-void SelectControlPoints(const Mat3X &X_world, 
-                         Mat *X_centered, 
-                         Mat34 *X_control_points) {
+static void SelectControlPoints(const Mat3X &X_world,
+                                Mat *X_centered,
+                                Mat34 *X_control_points) {
   size_t num_points = X_world.cols();
 
   // The first virtual control point, C0, is the centroid.
@@ -377,9 +378,9 @@ void SelectControlPoints(const Mat3X &X_world,
 }
 
 // Computes the barycentric coordinates for all real points
-void ComputeBarycentricCoordinates(const Mat3X &X_world_centered, 
-                                   const Mat34 &X_control_points,
-                                   Mat4X *alphas) {
+static void ComputeBarycentricCoordinates(const Mat3X &X_world_centered,
+                                          const Mat34 &X_control_points,
+                                          Mat4X *alphas) {
   size_t num_points = X_world_centered.cols();
   Mat3 C2 ;
   for (size_t c = 1; c < 4; c++) {
@@ -398,7 +399,7 @@ void ComputeBarycentricCoordinates(const Mat3X &X_world_centered,
 }
 
 // Estimates the coordinates of all real points in the camera coordinate frame
-void ComputePointsCoordinatesInCameraFrame(
+static void ComputePointsCoordinatesInCameraFrame(
     const Mat4X &alphas, 
     const Vec4 &betas,
     const Eigen::Matrix<double, 12, 12> &U,
@@ -435,8 +436,9 @@ void ComputePointsCoordinatesInCameraFrame(
 }
 
 bool EuclideanResectionEPnP(const Mat2X &x_camera,
-                            const Mat3X &X_world, 
-                            Mat3 *R, Vec3 *t) {
+                            const Mat3X &X_world,
+                            Mat3 *R, Vec3 *t,
+                            double success_threshold) {
   CHECK(x_camera.cols() == X_world.cols());
   CHECK(x_camera.cols() > 3);
   size_t num_points = X_world.cols();
@@ -535,7 +537,21 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
   vector<Vec3> ts(3);
   Vec rmse(3);
 
-  // TODO(julien): Document where the "1e-3" magical constant comes from below.
+  // At one point this threshold was 1e-3, and caused no end of problems in
+  // Blender by causing frames to not resect when they would have worked fine.
+  // When the resect failed, the projective followup is run leading to worse
+  // results, and often the dreaded "flipping" where objects get flipped
+  // between frames. Instead, disable the check for now, always succeeding. The
+  // ultimate check is always reprojection error anyway.
+  //
+  // TODO(keir): Decide if setting this to infinity, effectively disabling the
+  // check, is the right approach. So far this seems the case.
+  //
+  // TODO(sergey): Made it an option for now, in some cases it makes sense to
+  // still fallback to reprojection solution (see bug [#32765] from Blender bug tracker)
+
+  // double kSuccessThreshold = std::numeric_limits<double>::max();
+  double kSuccessThreshold = success_threshold;
 
   // Find the first possible solution for R, t corresponding to:
   // Betas          = [b00 b01 b11 b02 b12 b22 b03 b13 b23 b33]
@@ -548,7 +564,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
   Eigen::JacobiSVD<Mat> svd_of_l4(l_6x4, 
                                   Eigen::ComputeFullU | Eigen::ComputeFullV);
   Vec4 b4 = svd_of_l4.solve(rho);
-  if ((l_6x4 * b4).isApprox(rho, 1e-3)) {
+  if ((l_6x4 * b4).isApprox(rho, kSuccessThreshold)) {
     if (b4(0) < 0) {
       b4 = -b4;
     } 
@@ -574,7 +590,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
   Vec3 b3 = svdOfL3.solve(rho);
   VLOG(2) << " rho = " << rho;
   VLOG(2) << " l_6x3 * b3 = " << l_6x3 * b3;
-  if ((l_6x3 * b3).isApprox(rho, 1e-3)) {
+  if ((l_6x3 * b3).isApprox(rho, kSuccessThreshold)) {
     if (b3(0) < 0) {
       betas(0) = std::sqrt(-b3(0));
       betas(1) = (b3(2) < 0) ? std::sqrt(-b3(2)) : 0;
@@ -605,7 +621,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
   Eigen::JacobiSVD<Mat> svdOfL5(l_6x5, 
                                 Eigen::ComputeFullU | Eigen::ComputeFullV);
   Vec5 b5 = svdOfL5.solve(rho);
-  if ((l_6x5 * b5).isApprox(rho, 1e-3)) {
+  if ((l_6x5 * b5).isApprox(rho, kSuccessThreshold)) {
     if (b5(0) < 0) {
       betas(0) = std::sqrt(-b5(0));
       if (b5(2) < 0) {