hausdorff distance from a single triangle

Former-commit-id: 87256d2102b390843dc010fc85ee515d9135e97e

include/igl/copyleft/cgal/hausdorff.cpp

@@ -0,0 +1,39 @@
+#include "hausdorff.h"
+#include "../../hausdorff.h"
+#include <functional>
+
+template <
+  typename DerivedV,
+  typename Kernel,
+  typename Scalar>
+IGL_INLINE void igl::copyleft::cgal::hausdorff(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & treeB,
+  const std::vector<CGAL::Triangle_3<Kernel> > & /*TB*/,
+  Scalar & l,
+  Scalar & u)
+{
+  // Not sure why using `auto` here doesn't work with the `hausdorff` function
+  // parameter but explicitly naming the type does...
+  const std::function<double(const double &,const double &,const double &)> 
+    dist_to_B = [&treeB](
+    const double & x, const double & y, const double & z)->double
+  {
+    CGAL::Point_3<Kernel> query(x,y,z);
+    typename CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    >::Point_and_primitive_id pp = treeB.closest_point_and_primitive(query);
+    return std::sqrt((query-pp.first).squared_length());
+  };
+  return igl::hausdorff(V,dist_to_B,l,u);
+}

include/igl/copyleft/cgal/hausdorff.h

@@ -0,0 +1,62 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COPYLEFT_CGAL_HAUSDORFF_H
+#define IGL_COPYLEFT_CGAL_HAUSDORFF_H
+#include "../../igl_inline.h"
+
+#include <Eigen/Dense>
+#include "CGAL_includes.hpp"
+#include <vector>
+
+namespace igl 
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Compute lower and upper bounds (l,u) on the Hausdorff distance between a triangle
+      // (V) and a pointset (e.g., mesh, triangle soup) given by a distance function
+      // handle (dist_to_B).
+      //
+      // Inputs:
+      //   V   3 by 3 list of corner positions so that V.row(i) is the position of the
+      //     ith corner
+      //   treeB  CGAL's AABB tree containing triangle soup (VB,FB)
+      //   TB  list of CGAL triangles in order of FB (for determining which was found
+      //     in computation)
+      // Outputs:
+      //   l  lower bound on Hausdorff distance 
+      //   u  upper bound on Hausdorff distance
+      //
+      template <
+        typename DerivedV,
+        typename Kernel,
+        typename Scalar>
+      IGL_INLINE void hausdorff(
+        const Eigen::MatrixBase<DerivedV>& V,
+        const CGAL::AABB_tree<
+          CGAL::AABB_traits<Kernel, 
+            CGAL::AABB_triangle_primitive<Kernel, 
+              typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+            >
+          >
+        > & treeB,
+        const std::vector<CGAL::Triangle_3<Kernel> > & TB,
+        Scalar & l,
+        Scalar & u);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hausdorff.cpp"
+#endif
+
+#endif
+
+

include/igl/hausdorff.cpp

@@ -36,6 +36,54 @@ IGL_INLINE void igl::hausdorff(
   d = sqrt(std::max(dba,dab));
 }
 
+template <
+  typename DerivedV,
+  typename Scalar>
+IGL_INLINE void igl::hausdorff(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const std::function<Scalar(const Scalar &,const Scalar &, const Scalar &)> & dist_to_B,
+  Scalar & l,
+  Scalar & u)
+{
+  // e  3-long vector of opposite edge lengths
+  Eigen::Matrix<typename DerivedV::Scalar,1,3> e;
+  // Maximum edge length
+  Scalar e_max = 0;
+  for(int i=0;i<3;i++)
+  {
+    e(i) = (V.row((i+1)%3)-V.row((i+2)%3)).norm();
+    e_max = std::max(e_max,e(i));
+  }
+  // Semiperimeter
+  const Scalar s = (e(0)+e(1)+e(2))*0.5;
+  // Area
+  const Scalar A = sqrt(s*(s-e(0))*(s-e(1))*(s-e(2)));
+  // Circumradius
+  const Scalar R = e(0)*e(1)*e(2)/(4.*A);
+  // inradius
+  const Scalar r = A/s;
+  // Initialize lower bound to ∞
+  l = std::numeric_limits<Scalar>::infinity();
+  // d  3-long vector of distance from each corner to B
+  Eigen::Matrix<typename DerivedV::Scalar,1,3> d;
+  Scalar u1 = std::numeric_limits<Scalar>::infinity();
+  Scalar u2 = 0;
+  for(int i=0;i<3;i++)
+  {
+    d(i) = dist_to_B(V(i,0),V(i,1),V(i,2));
+    // Lower bound is simply the max over vertex distances
+    l = std::max(d(i),l);
+    // u1 is the minimum of corner distances + maximum adjacent edge 
+    u1 = std::min(u1,d(i) + std::max(e((i+1)%3),e((i+2)%3)));
+    // u2 first takes the maximum over corner distances
+    u2 = std::max(u2,d(i));
+  }
+  // u2 is the distance from the circumcenter/midpoint of obtuse edge plus the
+  // largest corner distance 
+  u2 += (s-r>2.*R ? R : 0.5*e_max);
+  u = std::min(u1,u2);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 template void igl::hausdorff<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>, double>(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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double&);
 #endif

include/igl/hausdorff.h

@@ -10,6 +10,7 @@
 #include "igl_inline.h"
 
 #include <Eigen/Dense>
+#include <functional>
 
 namespace igl 
 {
@@ -50,6 +51,28 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedVB> & VB, 
     const Eigen::PlainObjectBase<DerivedFB> & FB,
     Scalar & d);
+  // Compute lower and upper bounds (l,u) on the Hausdorff distance between a triangle
+  // (V) and a pointset (e.g., mesh, triangle soup) given by a distance function
+  // handle (dist_to_B).
+  //
+  // Inputs:
+  //   V   3 by 3 list of corner positions so that V.row(i) is the position of the
+  //     ith corner
+  //   dist_to_B  function taking the x,y,z coordinate of a query position and
+  //     outputing the closest-point distance to some point-set B
+  // Outputs:
+  //   l  lower bound on Hausdorff distance 
+  //   u  upper bound on Hausdorff distance
+  //
+  template <
+    typename DerivedV,
+    typename Scalar>
+  IGL_INLINE void hausdorff(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const std::function<
+      Scalar(const Scalar &,const Scalar &, const Scalar &)> & dist_to_B,
+    Scalar & l,
+    Scalar & u);
 }
 
 #ifndef IGL_STATIC_LIBRARY