Merge pull request #1630 from daviesthomas/master

Signing distances with Fast Winding Number

include/igl/fast_winding_number.cpp

@@ -434,21 +434,37 @@ IGL_INLINE void igl::fast_winding_number(
   assert(Q.cols() == 3 && "Q should be 3D");
   W.resize(Q.rows(),1);
   igl::parallel_for(Q.rows(),[&](int p)
-    {
+  {
     FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
-          Qp[0] = Q(p,0);
-          Qp[1] = Q(p,1);
-          Qp[2] = Q(p,2);
-      W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(
-        Qp,
-        accuracy_scale) 
-        / (4.0*igl::PI);
-    },1000);
+    Qp[0] = Q(p,0);
+    Qp[1] = Q(p,1);
+    Qp[2] = Q(p,2);
+    W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+  },1000);
+}
+
+template <typename Derivedp>
+IGL_INLINE typename Derivedp::Scalar igl::fast_winding_number(
+  const FastWindingNumberBVH & fwn_bvh,
+  const float accuracy_scale,
+  const Eigen::MatrixBase<Derivedp> & p)
+{
+  assert(p.cols() == 3 && "p should be 3D");
+
+  FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
+  Qp[0] = p(0,0);
+  Qp[1] = p(0,1);
+  Qp[2] = p(0,2);
+
+  typename Derivedp::Scalar w = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+
+  return w;
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&);
 // generated by autoexplicit.sh
 template void igl::fast_winding_number<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<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
@@ -457,4 +473,10 @@ template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
 template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, igl::FastWindingNumberBVH&);
 template void igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 template void igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, igl::FastWindingNumberBVH&);
+
+// tom did this manually. Unsure how to generate otherwise... sorry.
+template Eigen::Matrix<float, 1, 3, 1, 1, 3>::Scalar igl::fast_winding_number<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&);
+template void igl::fast_winding_number<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, int, igl::FastWindingNumberBVH&);
+template void igl::fast_winding_number<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int, igl::FastWindingNumberBVH&);
+template Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const>::Scalar igl::fast_winding_number<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const> > const&);
 #endif

include/igl/fast_winding_number.h

@@ -206,6 +206,19 @@ namespace igl
     const float accuracy_scale,
     const Eigen::MatrixBase<DerivedQ> & Q,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // After precomputation, compute winding number at a a single point
+  //
+  // Inputs:
+  //   fwn_bvh  Precomputed bounding volume hierarchy
+  //   accuracy_scale  parameter controlling accuracy (e.g., 2)
+  //   p single position
+  // Outputs:
+  //   w  winding number of this point
+  template <typename Derivedp>
+  IGL_INLINE typename Derivedp::Scalar fast_winding_number(
+    const FastWindingNumberBVH & fwn_bvh,
+    const float accuracy_scale,
+    const Eigen::MatrixBase<Derivedp> & p);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "fast_winding_number.cpp"

include/igl/signed_distance.cpp

@@ -13,6 +13,7 @@
 #include "per_vertex_normals.h"
 #include "point_mesh_squared_distance.h"
 #include "pseudonormal_test.h"
+#include "fast_winding_number.h"
 
 
 template <
@@ -37,10 +38,17 @@ IGL_INLINE void igl::signed_distance(
 {
   using namespace Eigen;
   using namespace std;
+
   const int dim = V.cols();
+
   assert((V.cols() == 3||V.cols() == 2) && "V should have 3d or 2d positions");
   assert((P.cols() == 3||P.cols() == 2) && "P should have 3d or 2d positions");
   assert(V.cols() == P.cols() && "V should have same dimension as P");
+
+  if (sign_type == SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER){
+    assert(V.cols() == 3 && "V should be 3D for fast winding number");
+  }
+
   // Only unsigned distance is supported for non-triangles
   if(sign_type != SIGNED_DISTANCE_TYPE_UNSIGNED)
   {
@@ -67,6 +75,9 @@ IGL_INLINE void igl::signed_distance(
   Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> E;
   Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> EMAP;
   WindingNumberAABB<RowVector3S,DerivedV,DerivedF> hier3;
+  igl::FastWindingNumberBVH fwn_bvh;
+  Eigen::VectorXf W;
+
   switch(sign_type)
   {
     default:
@@ -88,6 +99,10 @@ IGL_INLINE void igl::signed_distance(
           break;
       }
       break;
+    case SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER:
+      // assert above ensures dim == 3
+      igl::fast_winding_number(V.template cast<float>().eval(), F, 2, fwn_bvh);
+
     case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
       switch(dim)
       {
@@ -187,6 +202,14 @@ IGL_INLINE void igl::signed_distance(
           }
           break;
         }
+        case SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER:
+        {
+          //assert above ensured 3D
+          Scalar w = fast_winding_number(fwn_bvh, 2, q3.template cast<float>().eval());         
+          s = 1.-2.*std::abs(w);  
+
+          break;
+        }
         case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
         {
           RowVector3S n3;
@@ -439,7 +462,6 @@ IGL_INLINE void igl::signed_distance_winding_number(
   using namespace std;
   typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
   sqrd = tree.squared_distance(V,F,RowVector2S(q),i,(RowVector2S&)c);
-  Scalar w;
   // TODO: using .data() like this is very dangerous... This is assuming
   // colmajor order
   assert(!V.derived().IsRowMajor);
@@ -447,6 +469,55 @@ IGL_INLINE void igl::signed_distance_winding_number(
   s = 1.-2.*winding_number(V,F,q);
 }
 
+//Multi point by parrallel for on single point
+template <
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedS>
+IGL_INLINE void igl::signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh,
+    Eigen::PlainObjectBase<DerivedS> & S)
+  {
+    typedef Eigen::Matrix<typename DerivedV::Scalar,1,3> RowVector3S;
+    S.resize(P.rows(),1);
+    int min_parallel = 10000; 
+    parallel_for(P.rows(), [&](const int p)
+    {
+      RowVector3S q;
+      q.head(P.row(p).size()) = P.row(p);
+      // get sdf for single point, update result matrix
+      S(p) = signed_distance_fast_winding_number(q, V, F, tree,fwn_bvh);
+    }
+    ,min_parallel);  
+  }
+
+//Single Point
+template <
+  typename Derivedq,
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE typename DerivedV::Scalar igl::signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<Derivedq> & q,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh)
+  {
+    typedef typename DerivedV::Scalar Scalar;
+    Scalar s,sqrd;
+    Eigen::Matrix<Scalar,1,3> c;
+    int i = -1;
+    sqrd = tree.squared_distance(V,F,q,i,c);
+    Scalar w = fast_winding_number(fwn_bvh,2,q.template cast<float>());
+    //0.5 is on surface
+    return sqrt(sqrd)*(1.-2.*std::abs(w));
+  }
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
@@ -463,4 +534,5 @@ template void igl::signed_distance_pseudonormal<Eigen::Matrix<double, -1, -1, 0,
 template void igl::signed_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, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::SignedDistanceType, 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, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::signed_distance_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, double&, double&, int&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::signed_distance_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&);
+template void igl::signed_distance_fast_winding_number<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, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, igl::FastWindingNumberBVH const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/signed_distance.h

@@ -11,6 +11,7 @@
 #include "igl_inline.h"
 #include "AABB.h"
 #include "WindingNumberAABB.h"
+#include "fast_winding_number.h"
 #include <Eigen/Core>
 #include <vector>
 namespace igl
@@ -18,11 +19,14 @@ namespace igl
   enum SignedDistanceType
   {
     // Use fast pseudo-normal test [Bærentzen & Aanæs 2005]
-    SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
-    SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
-    SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
-    SIGNED_DISTANCE_TYPE_UNSIGNED       = 3,
-    NUM_SIGNED_DISTANCE_TYPE            = 4
+    SIGNED_DISTANCE_TYPE_PSEUDONORMAL         = 0,
+    // Use winding number [Jacobson, Kavan Sorking-Hornug 2013]
+    SIGNED_DISTANCE_TYPE_WINDING_NUMBER       = 1,
+    SIGNED_DISTANCE_TYPE_DEFAULT              = 2,
+    SIGNED_DISTANCE_TYPE_UNSIGNED             = 3,
+    // Use Fast winding number [Barill, Dickson, Schmidt, Levin, Jacobson 2018]
+    SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER  = 4,
+    NUM_SIGNED_DISTANCE_TYPE                  = 5
   };
   // Computes signed distance to a mesh
   //
@@ -62,7 +66,22 @@ namespace igl
     Eigen::PlainObjectBase<DerivedI> & I,
     Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedN> & N);
-  // Default bounds
+  // Computes signed distance to a mesh, with default bounds
+  //
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by ss list of triangle indices, ss should be 3 unless sign_type ==
+  //     SIGNED_DISTANCE_TYPE_UNSIGNED
+  //   sign_type  method for computing distance _sign_ S
+  //   lower_bound  lower bound of distances needed {std::numeric_limits::min}
+  //   upper_bound  lower bound of distances needed {std::numeric_limits::max}
+  // 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)
   template <
     typename DerivedP,
     typename DerivedV,
@@ -80,7 +99,7 @@ namespace igl
     Eigen::PlainObjectBase<DerivedI> & I,
     Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedN> & N);
-  // Computes signed distance to mesh
+  // Computes signed distance to mesh using pseudonormal with precomputed AABB tree and edge/vertice normals
   //
   // Inputs:
   //   tree  AABB acceleration tree (see AABB.h)
@@ -236,6 +255,63 @@ namespace igl
     Scalar & sqrd,
     int & i,
     Eigen::PlainObjectBase<Derivedc> & c);
+
+
+  // Calculates signed distance at query points P, using fast winding number
+  //   for sign.
+  //
+  // Usage:
+  //     VectorXd S;  
+  //     VectorXd V, P; //where V is mesh vertices, P are query points
+  //     VectorXi F;  
+  //     igl::FastWindingNumberBVH fwn_bvh;
+  //     igl::fast_winding_number(V.cast<float>(), F, 2, fwn_bvh);
+  //     igl::signed_distance_fast_winding_number(P,V,F,tree,fwn_bvh,S);
+  //
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of triangle indices
+  //   F  #F by 3 list of triangle normals 
+  //   tree  AABB acceleration tree (see AABB.h)
+  //   bvh fast winding precomputation (see Fast_Winding_Number.h)   
+  // Outputs:
+  //   S  #P list of signed distances of each point in P
+  template <
+    typename DerivedP,
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedS>
+  IGL_INLINE void signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh,
+    Eigen::PlainObjectBase<DerivedS> & S
+  );
+
+  // Calculates signed distance at query point q, using fast winding number
+  //   for sign.
+  //
+  // Inputs:
+  //   tree  AABB acceleration tree (see AABB.h)
+  //   V  #V by 3 list of triangle indices
+  //   F  #F by 3 list of triangle normals 
+  //   bvh fast winding precomputation (see Fast_Winding_Number.h)   
+  //   q  1 by 3 list of query point positions
+  // Outputs:
+  //   S  #P list of signed distances of each point in P
+  template <
+    typename Derivedq,
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE typename DerivedV::Scalar signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<Derivedq> & q,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh
+  );
 }
 
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/signed_distance.cpp

@@ -15,6 +15,7 @@ TEST_CASE("signed_distance: single_tet", "[igl]")
     0,2,1,
     0,3,2,
     1,2,3;
+
   Eigen::MatrixXd P(1,3);
   P<<0.5,0.5,0.5;
   for(const igl::SignedDistanceType type :
@@ -22,14 +23,23 @@ TEST_CASE("signed_distance: single_tet", "[igl]")
       igl::SIGNED_DISTANCE_TYPE_PSEUDONORMAL  ,
       igl::SIGNED_DISTANCE_TYPE_WINDING_NUMBER,
       igl::SIGNED_DISTANCE_TYPE_DEFAULT       ,
-      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      })
+      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      ,
+      igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER 
+      })
   {
     Eigen::VectorXd S;
     Eigen::VectorXi I;
     Eigen::MatrixXd C,N;
     igl::signed_distance( P,V,F,type,S,I,C,N);
     Eigen::VectorXd Sexact (1,1);Sexact<<sqrt(1./12.);
-    test_common::assert_near(S,Sexact,1e-15);
+
+    if (type == igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER) {
+      // loosen tolerance on fwn. 
+      test_common::assert_near(S,Sexact,1e-7);
+    } else {
+      test_common::assert_near(S,Sexact,1e-15);
+    }
+    
   }
 }
 
@@ -52,7 +62,7 @@ TEST_CASE("signed_distance: single_triangle", "[igl]")
       igl::SIGNED_DISTANCE_TYPE_PSEUDONORMAL  ,
       igl::SIGNED_DISTANCE_TYPE_WINDING_NUMBER,
       igl::SIGNED_DISTANCE_TYPE_DEFAULT       ,
-      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      
+      igl::SIGNED_DISTANCE_TYPE_UNSIGNED
       })
   {
     Eigen::VectorXd S;
@@ -63,3 +73,4 @@ TEST_CASE("signed_distance: single_triangle", "[igl]")
     test_common::assert_near(S,Sexact,1e-15);
   }
 }
+

tutorial/704_SignedDistance/main.cpp

@@ -19,7 +19,10 @@
 
 Eigen::MatrixXd V;
 Eigen::MatrixXi T,F;
+
 igl::AABB<Eigen::MatrixXd,3> tree;
+igl::FastWindingNumberBVH fwn_bvh;
+
 Eigen::MatrixXd FN,VN,EN;
 Eigen::MatrixXi E;
 Eigen::VectorXi EMAP;
@@ -28,6 +31,8 @@ double max_distance = 1;
 double slice_z = 0.5;
 bool overlay = false;
 
+bool useFastWindingNumber = false;
+
 void update_visualization(igl::opengl::glfw::Viewer & viewer)
 {
   using namespace Eigen;
@@ -72,12 +77,16 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
 
   // Compute signed distance
   VectorXd S_vis;
+
+  if (!useFastWindingNumber)
   {
     VectorXi I;
     MatrixXd N,C;
     // Bunny is a watertight mesh so use pseudonormal for signing
     signed_distance_pseudonormal(V_vis,V,F,tree,FN,VN,EN,EMAP,S_vis,I,C,N);
-  }
+  } else {
+    signed_distance_fast_winding_number(V_vis, V, F, tree, fwn_bvh, S_vis);
+  }    
 
   const auto & append_mesh = [&F_vis,&V_vis](
     const Eigen::MatrixXd & V,
@@ -114,6 +123,12 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
     case ',':
       slice_z = std::max(slice_z-0.01,0.01);
       break;
+    case '1':
+      useFastWindingNumber = true;
+      break;
+    case '2':
+      useFastWindingNumber = false;
+      break;
   }
   update_visualization(viewer);
   return true;
@@ -127,6 +142,7 @@ int main(int argc, char *argv[])
   cout<<"Usage:"<<endl;
   cout<<"[space]  toggle showing surface."<<endl;
   cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
+  cout<< "1/2 toggle between fast winding number (1) and pseudonormal (2) signing. \n";
   cout<<endl;
 
   // Load mesh: (V,T) tet-mesh of convex hull, F contains original surface
@@ -143,6 +159,9 @@ int main(int argc, char *argv[])
     max_distance = sqrt(sqrD.maxCoeff());
   }
 
+  // Fast winding and Pseudo normal depend on differnt AABB trees... We initialize both here.
+
+  // Pseudonormal setup...
   // Precompute signed distance AABB tree
   tree.init(V,F);
   // Precompute vertex,edge and face normals
@@ -152,6 +171,9 @@ int main(int argc, char *argv[])
   igl::per_edge_normals(
     V,F,igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
 
+  // fast winding number setup (just init fwn bvh)
+  igl::fast_winding_number(V, F, 2, fwn_bvh);
+
   // Plot the generated mesh
   igl::opengl::glfw::Viewer viewer;
   update_visualization(viewer);