revert hausdorff

Former-commit-id: 70b5f92dff0852d5adc8246fdf8b8840857ddf97

include/igl/hausdorff.cpp

@@ -7,116 +7,6 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "hausdorff.h"
 #include "point_mesh_squared_distance.h"
-#include "AABB.h"
-#include "upsample.h"
-#include "slice_mask.h"
-#include "remove_unreferenced.h"
-#include "EPS.h"
-
-template <
-  typename DerivedVA, 
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename Scalar>
-IGL_INLINE void igl::hausdorff(
-  const Eigen::PlainObjectBase<DerivedVA> & OVA, 
-  const Eigen::PlainObjectBase<DerivedFA> & OFA,
-  const Eigen::PlainObjectBase<DerivedVB> & VB, 
-  const Eigen::PlainObjectBase<DerivedFB> & FB,
-  const Scalar eps,
-  Scalar & lower,
-  Scalar & upper)
-{
-  using namespace Eigen;
-  using namespace igl;
-  using namespace std;
-  assert(OVA.cols() == 3 && "VA should contain 3d points");
-  assert(OFA.cols() == 3 && "FA should contain triangles");
-  assert(VB.cols() == 3 && "VB should contain 3d points");
-  assert(FB.cols() == 3 && "FB should contain triangles");
-  assert(eps>0 && "eps should be a strictly positive number");
-  const auto inf = std::numeric_limits<Scalar>::infinity();
-  DerivedVA VA = OVA;
-  DerivedFA FA = OFA;
-  // Precompute acceleration data structure for B
-  AABB<DerivedVB,3> tree;
-  tree.init(VB,FB);
-  lower = 0;
-  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
-  VectorXS DV;
-  Eigen::Matrix<int,Eigen::Dynamic,1> DI;
-  Eigen::Matrix<Scalar,Eigen::Dynamic,3> DC;
-  while(true)
-  {
-    // Compute distance to each vertex in VA
-    tree.squared_distance(VB,FB,VA,DV,DI,DC);
-    DV = DV.array().sqrt();
-    // Compute upper bound for each face in FA
-    VectorXS U(FA.rows(),1);
-    for(int i = 0;i<FA.rows();i++)
-    {
-      // Compute edge lengths
-      Eigen::Matrix<Scalar,3,1> e;
-      for(int c = 0;c<3;c++)
-      {
-        e(c) = (VA.row(FA(i,(c+1)%3)) - VA.row(FA(i,(c+2)%3))).norm();
-      }
-      // Semiperimeter
-      Scalar s = e.array().sum()/2.0;
-      // Area
-      Scalar A = sqrt(s*(s-e(0))*(s-e(1))*(s-e(2)));
-      // Circumradius
-      Scalar R = e(0)*e(1)*e(2)/(4.0*A);
-      // Inradius
-      Scalar r = A/s;
-      // Barycenter
-      Eigen::Matrix<typename DerivedVA::Scalar,1,3> B = 
-        (VA.row(FA(i,0))+ VA.row(FA(i,1))+ VA.row(FA(i,2)))/3.0;
-      // These still sometimes contribute
-      Scalar u1 = inf;
-      Scalar u2 = 0;
-      // u3 is a _huge_ improvement over u1, u2
-      Scalar u3 = 0;
-      // u4 _is contributing_ but it's rather expensive
-      Scalar u4 = inf;
-      for(int c = 0;c<3;c++)
-      {
-        const Scalar dic = DV(FA(i,c));
-        lower = std::max(lower,dic);
-        u1 = std::min(u1,dic+max(e((c+1)%3), e((c+2)%3)));
-        u2 = std::max(u2,dic);
-        u3 = std::max(u3,dic);
-        u3 = std::max(u3, (B-DC.row(FA(i,c))).norm());
-        {
-          Scalar u4c = 0;
-          for(int o = 0;o<3;o++)
-          {
-            u4c = 
-              std::max(u4c,c==o?dic:(VA.row(FA(i,c))-DC.row(FA(i,o))).norm());
-          }
-          u4 = std::min(u4,u4c);
-        }
-      }
-      u2 += ( s-r > 2.*R ? R : e.maxCoeff()/2.0 );
-      U(i) = std::min(std::min(std::min(u1,u2),u3),u4);
-    }
-    upper = U.maxCoeff();
-    //cout<<std::setprecision(17) << FA.rows()<<"\t"<<lower<< "\t"<<upper<<endl;
-    if(upper-lower <= eps)
-    {
-      return;
-    }
-    // Remove faces with too small upper bounds to impact result
-    slice_mask(DerivedFA(FA),(U.array()>lower).eval(),1,FA);
-    {
-      Eigen::Matrix<typename DerivedFA::Index,Eigen::Dynamic,1> _;
-      remove_unreferenced(DerivedVA(VA),DerivedFA(FA),VA,FA,_);
-    }
-    // Upsample all (remaining) faces in A
-    upsample(DerivedVA(VA),DerivedFA(FA),VA,FA);
-  }
-}
 
 template <
   typename DerivedVA, 
@@ -132,14 +22,18 @@ IGL_INLINE void igl::hausdorff(
   Scalar & d)
 {
   using namespace Eigen;
-  auto X = VA.colwise().maxCoeff().array().max(VB.colwise().maxCoeff().array());
-  auto N = VA.colwise().minCoeff().array().min(VB.colwise().minCoeff().array());
-  const double bbd = (X-N).matrix().norm();
-  const double eps = igl::EPS<Scalar>()*bbd;
-  double lab,uab,lba,uba;
-  igl::hausdorff(VA,FA,VB,FB,eps,lab,uab);
-  igl::hausdorff(VB,FB,VA,FA,eps,lba,uba);
-  d = std::max(lab,lba);
+  assert(VA.cols() == 3 && "VA should contain 3d points");
+  assert(FA.cols() == 3 && "FA should contain triangles");
+  assert(VB.cols() == 3 && "VB should contain 3d points");
+  assert(FB.cols() == 3 && "FB should contain triangles");
+  Matrix<Scalar,Dynamic,1> sqr_DBA,sqr_DAB;
+  Matrix<typename DerivedVA::Index,Dynamic,1> I;
+  Matrix<typename DerivedVA::Scalar,Dynamic,3> C;
+  point_mesh_squared_distance(VB,VA,FA,sqr_DBA,I,C);
+  point_mesh_squared_distance(VA,VB,FB,sqr_DAB,I,C);
+  const Scalar dba = sqr_DBA.maxCoeff();
+  const Scalar dab = sqr_DAB.maxCoeff();
+  d = sqrt(std::max(dba,dab));
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/hausdorff.h

@@ -13,51 +13,20 @@
 
 namespace igl 
 {
-  // HAUSDORFF compute lower and upper bounds on the **directed** (asymmetric)
-  // Hausdorff distance from mesh (VA,FA) to mesh (VB,FB) within a given error
-  // bound epsilon:
-  //
-  // l(A,B) ≤ max min d(a,b) ≤ u(A,B)
-  //          a∈A b∈B
-  // u(A,B)-l(A,B) ≤ epsilon
-  //
-  // This is mostly an implementation of "Fast and accurate Hausdorff distance
-  // calculation between meshes" [Guthe et al. 2005], with some ideas from
-  // "Interactive Hausdorff Distance Computation for General Polygonal Models"
-  // [Tang et al. 2009].
-  // 
-  // Inputs:
-  //   VA  #VA by 3 list of vertex positions
-  //   FA  #FA by 3 list of face indices into VA
-  //   VB  #VB by 3 list of vertex positions
-  //   FB  #FB by 3 list of face indices into VB
-  //   eps  desired bound on error
-  // Outputs:
-  //   lower  lower bound on directed Hausdorff distance
-  //   upper  upper bound on directed Hausdorff distance
-  template <
-    typename DerivedVA, 
-    typename DerivedFA,
-    typename DerivedVB,
-    typename DerivedFB,
-    typename Scalar>
-  IGL_INLINE void hausdorff(
-    const Eigen::PlainObjectBase<DerivedVA> & OVA, 
-    const Eigen::PlainObjectBase<DerivedFA> & OFA,
-    const Eigen::PlainObjectBase<DerivedVB> & VB, 
-    const Eigen::PlainObjectBase<DerivedFB> & FB,
-    const Scalar eps,
-    Scalar & lower,
-    Scalar & upper);
-  // HAUSDORFF compute the **symmetric** Hausdorff distance between mesh
-  // (VA,FA) and mesh (VB,FB). That is:
+  // HAUSDORFF compute the Hausdorff distance between mesh (VA,FA) and mesh
+  // (VB,FB). This is the 
   //
   // d(A,B) = max ( max min d(a,b) , max min d(b,a) )
   //                a∈A b∈B          b∈B a∈A
   //
-  // Note: this is an iterative method that will compute the distance within
-  // double precision epsilon times the bounding box diagonal of the combined
-  // meshes A and B.
+  // Known issue: This is only computing max(min(va,B),min(vb,A)). This is
+  // better than max(min(va,Vb),min(vb,Va)). This (at least) is missing
+  // "edge-edge" cases like the distance between the two different
+  // triangulations of a non-planar quad in 3D. Even simpler, consider the
+  // Hausdorff distance between the non-convex, block letter V polygon (with 7
+  // vertices) in 2D and its convex hull. The Hausdorff distance is defined by
+  // the midpoint in the middle of the segment across the concavity and some
+  // non-vertex point _on the edge_ of the V.
   //
   // Inputs:
   //   VA  #VA by 3 list of vertex positions
@@ -66,6 +35,8 @@ namespace igl
   //   FB  #FB by 3 list of face indices into VB
   // Outputs:
   //   d  hausdorff distance
+  //   //pair  2 by 3 list of "determiner points" so that pair(1,:) is from A
+  //   //  and pair(2,:) is from B
   //
   template <
     typename DerivedVA,