pseudonormal test and signed distance in 2D

Former-commit-id: 1c9af2807cd3b737cda9169a428059e82535cf1e

include/igl/pseudonormal_test.cpp

@@ -62,3 +62,44 @@ IGL_INLINE void igl::pseudonormal_test(
   }
   s = (qc.dot(n) >= 0 ? 1. : -1.);
 }
+
+IGL_INLINE void igl::pseudonormal_test(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & E,
+  const Eigen::MatrixXd & EN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::RowVector2d & q,
+  const int e,
+  const Eigen::RowVector2d & c,
+  double & s,
+  Eigen::RowVector2d & n)
+{
+  using namespace Eigen;
+  const auto & qc = q-c;
+  const double len = (V.row(E(e,1))-V.row(E(e,0))).norm();
+  // barycentric coordinates
+  RowVector2d b((c-V.row(E(e,1))).norm()/len,(c-V.row(E(e,0))).norm()/len);
+  // Determine which normal to use
+  const double epsilon = 1e-12;
+  const int type = (b.array()<=epsilon).cast<int>().sum();
+  switch(type)
+  {
+    case 1:
+      // Find vertex
+      for(int x = 0;x<2;x++)
+      {
+        if(b(x)>epsilon)
+        {
+          n = VN.row(E(e,x));
+          break;
+        }
+      }
+      break;
+    default:
+      assert(false && "all barycentric coords zero.");
+    case 0:
+      n = EN.row(e);
+      break;
+  }
+  s = (qc.dot(n) >= 0 ? 1. : -1.);
+}

include/igl/pseudonormal_test.h

@@ -41,6 +41,16 @@ namespace igl
     const Eigen::RowVector3d & c,
     double & s,
     Eigen::RowVector3d & n);
+  IGL_INLINE void pseudonormal_test(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & E,
+    const Eigen::MatrixXd & EN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::RowVector2d & q,
+    const int i,
+    const Eigen::RowVector2d & c,
+    double & s,
+    Eigen::RowVector2d & n);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "pseudonormal_test.cpp"

include/igl/signed_distance.cpp

@@ -26,22 +26,34 @@ IGL_INLINE void igl::signed_distance(
 {
   using namespace Eigen;
   using namespace std;
-  assert(V.cols() == 3 && "V should have 3d positions");
-  assert(P.cols() == 3 && "P should have 3d positions");
+  const int dim = V.cols();
+  assert((V.cols() == 3||V.cols() == 2) && "V should have 3d or 2d positions");
+  assert((P.cols() == 3||V.cols() == 2) && "P should have 3d or 2d positions");
+  assert(V.cols() == P.cols() && "V should have same dimension as P");
   // Only unsigned distance is supported for non-triangles
   if(sign_type != SIGNED_DISTANCE_TYPE_UNSIGNED)
   {
-    assert(F.cols() == 3 && "F should have triangles");
+    assert(F.cols() == dim && "F should have co-dimension 0 simplices");
   }
 
   // Prepare distance computation
-  AABB<MatrixXd,3> tree;
-  tree.init(V,F);
+  AABB<MatrixXd,3> tree3;
+  AABB<MatrixXd,2> tree2;
+  switch(dim)
+  {
+    default:
+    case 3:
+      tree3.init(V,F);
+      break;
+    case 2:
+      tree2.init(V,F);
+      break;
+  }
 
   Eigen::MatrixXd FN,VN,EN;
   Eigen::MatrixXi E;
   Eigen::VectorXi EMAP;
-  WindingNumberAABB<Eigen::Vector3d> hier;
+  WindingNumberAABB<Eigen::Vector3d> hier3;
   switch(sign_type)
   {
     default:
@@ -51,28 +63,63 @@ IGL_INLINE void igl::signed_distance(
       break;
     case SIGNED_DISTANCE_TYPE_DEFAULT:
     case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-      hier.set_mesh(V,F);
-      hier.grow();
+      switch(dim)
+      {
+        default:
+        case 3:
+          hier3.set_mesh(V,F);
+          hier3.grow();
+          break;
+        case 2:
+          // no precomp, no hierarchy
+          break;
+      }
       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,FN,EN,E,EMAP);
-      N.resize(P.rows(),3);
+      switch(dim)
+      {
+        default:
+        case 3:
+          // "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,FN,EN,E,EMAP);
+          break;
+        case 2:
+          FN.resize(F.rows(),2);
+          VN = MatrixXd::Zero(V.rows(),2);
+          for(int e = 0;e<F.rows();e++)
+          {
+            // rotate edge vector
+            FN(e,0) = -(V(F(e,1),1)-V(F(e,0),1));
+            FN(e,1) =  (V(F(e,1),0)-V(F(e,0),0));
+            FN.row(e).normalize();
+            // add to vertex normal
+            VN.row(F(e,1)) += FN.row(e);
+            VN.row(F(e,0)) += FN.row(e);
+          }
+          // normalize to average
+          VN.rowwise().normalize();
+          break;
+      }
+      N.resize(P.rows(),dim);
       break;
   }
 
   S.resize(P.rows(),1);
   I.resize(P.rows(),1);
-  C.resize(P.rows(),3);
+  C.resize(P.rows(),dim);
   for(int p = 0;p<P.rows();p++)
   {
-    const RowVector3d q = P.row(p);
+    const RowVectorXd q = P.row(p);
+    const RowVector3d q3 = q;
+    const RowVector2d q2 = q;
     double s,sqrd;
-    RowVector3d c;
+    RowVectorXd c;
+    RowVector3d c3;
+    RowVector2d c2;
     int i=-1;
     switch(sign_type)
     {
@@ -80,23 +127,32 @@ IGL_INLINE void igl::signed_distance(
         assert(false && "Unknown SignedDistanceType");
       case SIGNED_DISTANCE_TYPE_UNSIGNED:
         s = 1.;
-        sqrd = tree.squared_distance(V,F,q,i,c);
+        sqrd = dim==3?
+          tree3.squared_distance(V,F,q3,i,c3):
+          tree2.squared_distance(V,F,q2,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_DEFAULT:
       case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-        signed_distance_winding_number(tree,V,F,hier,q,s,sqrd,i,c);
+        dim==3 ? 
+          signed_distance_winding_number(tree3,V,F,hier3,q3,s,sqrd,i,c3):
+          signed_distance_winding_number(tree2,V,F,q2,s,sqrd,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
       {
-        Eigen::RowVector3d n;
-        signed_distance_pseudonormal(tree,V,F,FN,VN,EN,EMAP,q,s,sqrd,i,c,n);
+        RowVector3d n3;
+        RowVector2d n2;
+        dim==3 ?
+          signed_distance_pseudonormal(tree3,V,F,FN,VN,EN,EMAP,q3,s,sqrd,i,c3,n3):
+          signed_distance_pseudonormal(tree2,V,F,FN,VN,q2,s,sqrd,i,c2,n2);
+        Eigen::RowVectorXd n;
+        (dim==3 ? n = n3 : n = n3);
         N.row(p) = n;
         break;
       }
     }
     I(p) = i;
     S(p) = s*sqrt(sqrd);
-    C.row(p) = c;
+    C.row(p) = (dim==3 ? c=c3 : c=c2);
   }
 }
 
@@ -191,6 +247,25 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   pseudonormal_test(V,F,FN,VN,EN,EMAP,q,f,c,s,n);
 }
 
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const AABB<Eigen::MatrixXd,2> & tree,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::RowVector2d & q,
+  double & s,
+  double & sqrd,
+  int & f,
+  Eigen::RowVector2d & c,
+  Eigen::RowVector2d & n)
+{
+  using namespace Eigen;
+  using namespace std;
+  sqrd = tree.squared_distance(V,F,q,f,c);
+  pseudonormal_test(V,F,FN,VN,q,f,c,s,n);
+}
+
 IGL_INLINE double igl::signed_distance_winding_number(
   const AABB<Eigen::MatrixXd,3> & tree,
   const Eigen::MatrixXd & V,
@@ -210,12 +285,12 @@ IGL_INLINE void igl::signed_distance_winding_number(
   const AABB<Eigen::MatrixXd,3> & tree,
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
-  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-  const Eigen::RowVector3d & q,
+  const igl::WindingNumberAABB<Eigen::Matrix<double,3,1> > & hier,
+  const Eigen::Matrix<double,1,3> & q,
   double & s,
   double & sqrd,
   int & i,
-  Eigen::RowVector3d & c)
+  Eigen::Matrix<double,1,3> & c)
 {
   using namespace Eigen;
   using namespace std;
@@ -225,6 +300,26 @@ IGL_INLINE void igl::signed_distance_winding_number(
   s = 1.-2.*w;
 }
 
+IGL_INLINE void igl::signed_distance_winding_number(
+  const AABB<Eigen::MatrixXd,2> & tree,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::Matrix<double,1,2> & q,
+  double & s,
+  double & sqrd,
+  int & i,
+  Eigen::Matrix<double,1,2> & c)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  sqrd = tree.squared_distance(V,F,q,i,c);
+  double w;
+  winding_number_2(V.data(), V.rows(), F.data(), F.rows(), q.data(), 1, &w);
+  s = 1.-2.*w;
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // This template is necessary for the others to compile with clang
 // http://stackoverflow.com/questions/27748442/is-clangs-c11-support-reliable

include/igl/signed_distance.h

@@ -77,6 +77,19 @@ namespace igl
   //   i  closest primitive
   //   c  closest point
   //   n  normal
+  IGL_INLINE void signed_distance_pseudonormal(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const AABB<Eigen::MatrixXd,3> & tree,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    Eigen::VectorXd & S,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C,
+    Eigen::MatrixXd & N);
   IGL_INLINE void signed_distance_pseudonormal(
     const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & V,
@@ -92,18 +105,17 @@ namespace igl
     Eigen::RowVector3d & c,
     Eigen::RowVector3d & n);
   IGL_INLINE void signed_distance_pseudonormal(
-    const Eigen::MatrixXd & P,
+    const AABB<Eigen::MatrixXd,2> & tree,
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & FN,
     const Eigen::MatrixXd & VN,
-    const Eigen::MatrixXd & EN,
-    const Eigen::VectorXi & EMAP,
-    Eigen::VectorXd & S,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C,
-    Eigen::MatrixXd & N);
+    const Eigen::RowVector2d & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::RowVector2d & c,
+    Eigen::RowVector2d & n);
 
   // Inputs:
   //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)
@@ -124,12 +136,21 @@ namespace igl
     const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-    const Eigen::RowVector3d & q,
+    const igl::WindingNumberAABB<Eigen::Matrix<double,3,1> > & hier,
+    const Eigen::Matrix<double,1,3> & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::Matrix<double,1,3> & c);
+  IGL_INLINE void signed_distance_winding_number(
+    const AABB<Eigen::MatrixXd,2> & tree,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::Matrix<double,1,2> & q,
     double & s,
     double & sqrd,
     int & i,
-    Eigen::RowVector3d & c);
+    Eigen::Matrix<double,1,2> & c);
 }
 
 #ifndef IGL_STATIC_LIBRARY