signed distance overloads

Former-commit-id: 94842b3e0ae758588e6849621af93da5502bee51

include/igl/WindingNumberAABB.h

@@ -31,6 +31,7 @@ namespace igl
         NUM_SPLIT_METHODS = 3
       } split_method;
     public:
+      inline WindingNumberAABB(){}
       inline WindingNumberAABB(
         const Eigen::MatrixXd & V,
         const Eigen::MatrixXi & F);
@@ -38,6 +39,9 @@ namespace igl
         const WindingNumberTree<Point> & parent,
         const Eigen::MatrixXi & F);
       // Initialize some things
+      inline void set_mesh(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
       inline void init();
       inline bool inside(const Point & p) const;
       inline virtual void grow();
@@ -67,6 +71,15 @@ namespace igl
 #  define WindingNumberAABB_MIN_F 100
 #endif
 
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::set_mesh(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F)
+{
+  igl::WindingNumberTree<Point>::set_mesh(V,F);
+  init();
+}
+
 template <typename Point>
 inline void igl::WindingNumberAABB<Point>::init()
 {

include/igl/WindingNumberTree.h

@@ -12,6 +12,7 @@
 #include <Eigen/Dense>
 #include "WindingNumberMethod.h"
 
+static Eigen::MatrixXd dummyV;
 namespace igl
 {
   // Templates:
@@ -44,6 +45,7 @@ namespace igl
       // (Approximate) center (of mass)
       Point center;
     public:
+      inline WindingNumberTree():V(dummyV){}
       // For root
       inline WindingNumberTree(
         const Eigen::MatrixXd & V,
@@ -53,6 +55,9 @@ namespace igl
         const WindingNumberTree<Point> & parent,
         const Eigen::MatrixXi & F);
       inline virtual ~WindingNumberTree();
+      inline virtual void set_mesh(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
       // Set method
       inline void set_method( const WindingNumberMethod & m);
     public:
@@ -137,7 +142,6 @@ template <typename Point>
 std::map< std::pair<const igl::WindingNumberTree<Point>*,const igl::WindingNumberTree<Point>*>, double>
   igl::WindingNumberTree<Point>::cached;
 
-static Eigen::MatrixXd dummyV;
 template <typename Point>
 inline igl::WindingNumberTree<Point>::WindingNumberTree(
   const Eigen::MatrixXd & _V,
@@ -151,6 +155,14 @@ inline igl::WindingNumberTree<Point>::WindingNumberTree(
   cap(),
   radius(std::numeric_limits<double>::infinity()),
   center(0,0,0)
+{
+  set_mesh(V,F);
+}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::set_mesh(
+    const Eigen::MatrixXd & _V,
+    const Eigen::MatrixXi & _F)
 {
   // Remove any exactly duplicate vertices
   // Q: Can this ever increase the complexity of the boundary?

include/igl/cgal/signed_distance.cpp

@@ -6,6 +6,100 @@
 // 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/.
 #include "signed_distance.h"
+#include "../per_vertex_normals.h"
+#include "../per_edge_normals.h"
+#include "../per_face_normals.h"
+#include "point_mesh_squared_distance.h"
+
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance(
+  const Eigen::MatrixXd & P,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const SignedDistanceType sign_type,
+  Eigen::VectorXd & S,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C,
+  Eigen::MatrixXd & N)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(V.cols() == 3 && "V should have 3d positions");
+  assert(P.cols() == 3 && "P should have 3d positions");
+  assert(F.cols() == 3 && "F should have triangles");
+
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  // Prepare distance computation
+  Tree tree;
+  vector<Triangle_3 > T;
+  point_mesh_squared_distance_precompute(V,F,tree,T);
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  WindingNumberAABB<Eigen::Vector3d> hier;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      hier.set_mesh(V,F);
+      hier.grow();
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+      // [Bærentzen & Aanæs 2005]
+      per_face_normals(V,F,FN);
+      per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+      per_edge_normals(
+        V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
+      N.resize(P.rows(),3);
+      break;
+  }
+
+  S.resize(P.rows(),1);
+  I.resize(P.rows(),1);
+  C.resize(P.rows(),3);
+  for(int p = 0;p<P.rows();p++)
+  {
+    const Point_3 q(P(p,0),P(p,1),P(p,2));
+    double s,sqrd;
+    Point_and_primitive_id pp;
+    switch(sign_type)
+    {
+      default:
+        assert(false && "Unknown SignedDistanceType");
+      case SIGNED_DISTANCE_TYPE_DEFAULT:
+      case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+        signed_distance_winding_number(tree,hier,q,s,sqrd,pp);
+        break;
+      case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      {
+        Vector3d n(0,0,0);
+        signed_distance_pseudonormal(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
+        N.row(p) = n;
+        break;
+      }
+    }
+    I(p) = pp.second - T.begin();
+    S(p) = s*sqrt(sqrd);
+    C(p,0) = pp.first.x();
+    C(p,1) = pp.first.y();
+    C(p,2) = pp.first.z();
+  }
+}
+
+
 template <typename Kernel>
 IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   const CGAL::AABB_tree<
@@ -22,6 +116,44 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   const Eigen::MatrixXd & EN,
   const Eigen::VectorXi & EMAP,
   const typename Kernel::Point_3 & q)
+{
+  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;
+  typename Kernel::FT s,sqrd;
+  Eigen::Vector3d n(0,0,0);
+  signed_distance_pseudonormal<Kernel>(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
+  return s*sqrt(sqrd);
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::MatrixXd & EN,
+  const Eigen::VectorXi & EMAP,
+  const typename Kernel::Point_3 & q,
+  typename Kernel::FT & s,
+  typename Kernel::FT & sqrd,
+  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,
+   Eigen::Vector3d & n)
 {
   using namespace Eigen;
   using namespace std;
@@ -34,10 +166,10 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
   typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  pp = tree.closest_point_and_primitive(q);
   Point_3 & p = pp.first;
   const auto & qp = q-p;
-  const FT sqrd = qp.squared_length();
+  sqrd = qp.squared_length();
   Vector3d v(qp.x(),qp.y(),qp.z());
   const int f = pp.second - T.begin();
   const Triangle_3 & t = *pp.second;
@@ -49,7 +181,6 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   Vector3d b(area(1,2),area(2,0),area(0,1));
   b /= b.sum();
   // Determine which normal to use
-  Vector3d n;
   const double epsilon = 1e-12;
   const int type = (b.array()<=epsilon).cast<int>().sum();
   switch(type)
@@ -82,8 +213,7 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
       n = FN.row(f);
       break;
   }
-  const double s = (v.dot(n) >= 0 ? 1. : -1.);
-  return s*sqrt(sqrd);
+  s = (v.dot(n) >= 0 ? 1. : -1.);
 }
 
 template <typename Kernel>
@@ -97,6 +227,38 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
   > & tree,
   const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
   const typename Kernel::Point_3 & q)
+{
+  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;
+  typename Kernel::FT s,sqrd;
+  signed_distance_winding_number<Kernel>(tree,hier,q,s,sqrd,pp);
+  return s*sqrt(sqrd);
+}
+
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance_winding_number(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+  const typename Kernel::Point_3 & q,
+  typename Kernel::FT & s,
+  typename Kernel::FT & sqrd,
+  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)
 {
   typedef typename Kernel::FT FT;
   typedef typename Kernel::Point_3 Point_3;
@@ -109,14 +271,13 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
   using namespace Eigen;
   using namespace std;
   using namespace igl;
-  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  pp = tree.closest_point_and_primitive(q);
   Point_3 & p = pp.first;
   const auto & qp = q-p;
   const Vector3d eq(q.x(),q.y(),q.z()); 
-  const FT sqrd = qp.squared_length();
+  sqrd = qp.squared_length();
   const double w = hier.winding_number(eq);
-  const FT s = 1.-2.*w;
-  return s*sqrt(sqrd);
+  s = 1.-2.*w;
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/cgal/signed_distance.h

@@ -19,8 +19,38 @@ namespace igl
     SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
     SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
     SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
-    NUM_SIGNED_DISTANCE_TYPE            = 3
+    SIGNED_DISTANCE_TYPE_IGNORE         = 3,
+    NUM_SIGNED_DISTANCE_TYPE            = 4
   };
+  // Computes signed distance to a mesh
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  // Outputs:
+  //   S  #P list of smallest signed distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //   N  #P by 3 list of closest normals (only set if
+  //     sign_type=SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices are ignored.
+  template <typename Kernel>
+  IGL_INLINE void signed_distance(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const SignedDistanceType sign_type,
+    Eigen::VectorXd & S,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C,
+    Eigen::MatrixXd & N);
   // Computes signed distance to mesh
   //
   // Templates:
@@ -53,6 +83,37 @@ namespace igl
     const Eigen::MatrixXd & EN,
     const Eigen::VectorXi & EMAP,
     const typename Kernel::Point_3 & q);
+  // Outputs:
+  //   s  sign
+  //   sqrd  squared distance
+  //   pp  closest point and primitve
+  //   n  normal
+  template <typename Kernel>
+  IGL_INLINE void signed_distance_pseudonormal(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    const typename Kernel::Point_3 & q,
+    typename Kernel::FT & s,
+    typename Kernel::FT & sqrd,
+    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,
+   Eigen::Vector3d & n);
+
   // Inputs:
   //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
   //   hier  Winding number evaluation hierarchy
@@ -69,6 +130,29 @@ namespace igl
     > & tree,
     const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
     const typename Kernel::Point_3 & q);
+  // Outputs:
+  //   s  sign
+  //   sqrd  squared distance
+  //   pp  closest point and primitve
+  template <typename Kernel>
+  IGL_INLINE void signed_distance_winding_number(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+    const typename Kernel::Point_3 & q,
+    typename Kernel::FT & s,
+    typename Kernel::FT & sqrd,
+    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);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/signed_distance_isosurface.cpp

@@ -62,20 +62,33 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
   typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Eigen::MatrixXd FN,VN,EN;
-  Eigen::MatrixXi E;
-  Eigen::VectorXi EMAP;
-  // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
-  // [Bærentzen & Aanæs 2005]
-  per_face_normals(IV,IF,FN);
-  per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,VN);
-  per_edge_normals(IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
-  // Prepare distance computation
   Tree tree;
   vector<Triangle_3 > T;
   point_mesh_squared_distance_precompute(IV,IF,tree,T);
-  WindingNumberAABB<Eigen::Vector3d> hier(IV,IF);
-  hier.grow();
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  WindingNumberAABB<Eigen::Vector3d> hier;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      hier.set_mesh(IV,IF);
+      hier.grow();
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+      // [Bærentzen & Aanæs 2005]
+      per_face_normals(IV,IF,FN);
+      per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+      per_edge_normals(
+        IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
+      break;
+  }
+
   Tr tr;            // 3D-Delaunay triangulation
   C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
   // defining the surface

include/igl/per_edge_normals.cpp

@@ -8,6 +8,7 @@
 template <
   typename DerivedV, 
   typename DerivedF, 
+  typename DerivedFN, 
   typename DerivedN,
   typename DerivedE,
   typename DerivedEMAP>
@@ -15,9 +16,11 @@ IGL_INLINE void igl::per_edge_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   const PerEdgeNormalsWeightingType weighting,
+  const Eigen::PlainObjectBase<DerivedFN>& FN,
   Eigen::PlainObjectBase<DerivedN> & N,
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+
 {
   using namespace Eigen;
   using namespace std;
@@ -32,8 +35,6 @@ IGL_INLINE void igl::per_edge_normals(
   unique_simplices(allE,E,_,EMAP);
   // now sort(allE,2) == E(EMAP,:), that is, if EMAP(i) = j, then E.row(j) is
   // the undirected edge corresponding to the directed edge allE.row(i).
-  MatrixXd FN;
-  per_face_normals(V,F,FN);
 
   Eigen::VectorXd W(F.rows());
   switch(weighting)
@@ -60,7 +61,25 @@ IGL_INLINE void igl::per_edge_normals(
     }
   }
   N.rowwise().normalize();
-  
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedN,
+  typename DerivedE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::per_edge_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const PerEdgeNormalsWeightingType weighting,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  Eigen::PlainObjectBase<DerivedN> FN;
+  per_face_normals(V,F,FN);
+  return per_edge_normals(V,F,weighting,FN,N,E,EMAP);
 }
 
 template <

include/igl/per_edge_normals.h

@@ -31,6 +31,21 @@ namespace igl
   //   E  #E by 2 matrix of edge indices per row
   //   EMAP  #E by 1 matrix of indices from all edges to E
   //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedFN,
+    typename DerivedN,
+    typename DerivedE,
+    typename DerivedEMAP>
+  IGL_INLINE void per_edge_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerEdgeNormalsWeightingType weight,
+    const Eigen::PlainObjectBase<DerivedFN>& FN,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
   template <
     typename DerivedV, 
     typename DerivedF, 
@@ -63,4 +78,3 @@ namespace igl
 #endif
 
 #endif
-