fixed Hausdorff distance to converge to bounds on exact value

Former-commit-id: 7a74722e1e0325d4c237cfc1ea7a3c75f327c5ea

include/igl/AABB.cpp

@@ -955,6 +955,8 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;
+// generated by autoexplicit.sh
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;
 // generated by autoexplicit.sh
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;

include/igl/hausdorff.cpp

@@ -7,6 +7,11 @@
 // 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, 
@@ -15,25 +20,126 @@ template <
   typename DerivedFB,
   typename Scalar>
 IGL_INLINE void igl::hausdorff(
-  const Eigen::PlainObjectBase<DerivedVA> & VA, 
-  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVA> & OVA, 
+  const Eigen::PlainObjectBase<DerivedFA> & OFA,
   const Eigen::PlainObjectBase<DerivedVB> & VB, 
   const Eigen::PlainObjectBase<DerivedFB> & FB,
-  Scalar & d)
+  const Scalar eps,
+  Scalar & lower,
+  Scalar & upper)
 {
   using namespace Eigen;
+  using namespace igl;
+  using namespace std;
   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));
+  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, 
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename Scalar>
+IGL_INLINE void igl::hausdorff(
+  const Eigen::PlainObjectBase<DerivedVA> & VA, 
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVB> & VB, 
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  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);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/hausdorff.h

@@ -13,20 +13,51 @@
 
 namespace igl 
 {
-  // HAUSDORFF compute the Hausdorff distance between mesh (VA,FA) and mesh
-  // (VB,FB). This is the 
+  // 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:
   //
   // d(A,B) = max ( max min d(a,b) , max min d(b,a) )
   //                a∈A b∈B          b∈B a∈A
   //
-  // 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.
+  // 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.
   //
   // Inputs:
   //   VA  #VA by 3 list of vertex positions
@@ -35,8 +66,6 @@ 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, 

include/igl/remove_unreferenced.cpp

@@ -100,6 +100,8 @@ IGL_INLINE void igl::remove_unreferenced(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/slice.cpp

@@ -248,6 +248,8 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template Eigen::Matrix<double, -1, 1, 0, -1, 1>  igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::Matrix<double, -1, -1, 0, -1, -1>  igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::Matrix<double, -1, 1, 0, -1, 1>  igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);