better comments and signed distance prototype for multiple points

Former-commit-id: bfcf0a10e83f1055e4417de3d130210178e55534

include/igl/AABB.h

@@ -261,7 +261,7 @@ private:
 public:
       template <int SS>
       static
-      void barycentric_coordinates(
+      inline void barycentric_coordinates(
         const RowVectorDIMS & p, 
         const RowVectorDIMS & a, 
         const RowVectorDIMS & b, 

include/igl/in_element.h

@@ -44,10 +44,6 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedQ> & Q,
     const AABB<DerivedV,DIM> & aabb,
     Eigen::SparseMatrix<Scalar> & I);
-  //
-  // Example:
-  //   InElementAABB aabb;
-  //   aabb.init(V,Ele);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/point_mesh_squared_distance.cpp

@@ -29,4 +29,5 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+template void igl::point_mesh_squared_distance<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<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, 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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/pseudonormal_test.cpp

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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/.
+#include "pseudonormal_test.h"
+#include "AABB.h"
+#include <cassert>
+
+IGL_INLINE void igl::pseudonormal_test(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::MatrixXd & EN,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::RowVector3d & q,
+  const int f,
+  const Eigen::RowVector3d & c,
+  double & s,
+  Eigen::RowVector3d & n)
+{
+  using namespace Eigen;
+  const auto & qc = q-c;
+  RowVector3d b;
+  AABB<Eigen::MatrixXd,3>::barycentric_coordinates(
+    c,V.row(F(f,0)),V.row(F(f,1)),V.row(F(f,2)),b);
+  // Determine which normal to use
+  const double epsilon = 1e-12;
+  const int type = (b.array()<=epsilon).cast<int>().sum();
+  switch(type)
+  {
+    case 2:
+      // Find vertex
+      for(int x = 0;x<3;x++)
+      {
+        if(b(x)>epsilon)
+        {
+          n = VN.row(F(f,x));
+          break;
+        }
+      }
+      break;
+    case 1:
+      // Find edge
+      for(int x = 0;x<3;x++)
+      {
+        if(b(x)<=epsilon)
+        {
+          n = EN.row(EMAP(F.rows()*x+f));
+          break;
+        }
+      }
+      break;
+    default:
+      assert(false && "all barycentric coords zero.");
+    case 0:
+      n = FN.row(f);
+      break;
+  }
+  s = (qc.dot(n) >= 0 ? 1. : -1.);
+}

include/igl/pseudonormal_test.h

@@ -0,0 +1,48 @@
+// 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_PSEUDONORMAL_TEST_H
+#define IGL_PSEUDONORMAL_TEST_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Given a mesh (V,F), a query point q, and a point on (V,F) c, determine
+  // whether q is inside (V,F) --> s=-1 or outside (V,F) s=1, based on the
+  // sign of the dot product between (q-c) and n, where n is the normal _at c_,
+  // carefully chosen according to [Bærentzen & Aanæs 2005]
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  //   FN  #F by 3 list of triangle normals 
+  //   VN  #V by 3 list of vertex normals (ANGLE WEIGHTING)
+  //   EN  #E by 3 list of edge normals (UNIFORM WEIGHTING)
+  //   EMAP  #F*3 mapping edges in F to E
+  //   q  Query point
+  //   i  index into F to face to which c belongs
+  //   c  Point on (V,F)
+  // Outputs:
+  //   s  sign
+  //   n  normal
+  IGL_INLINE void pseudonormal_test(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::RowVector3d & q,
+    const int i,
+    const Eigen::RowVector3d & c,
+    double & s,
+    Eigen::RowVector3d & n);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "pseudonormal_test.cpp"
+#endif
+#endif

include/igl/signed_distance.cpp

@@ -6,11 +6,12 @@
 // 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 "get_seconds.h"
 #include "per_edge_normals.h"
 #include "per_face_normals.h"
-#include "get_seconds.h"
+#include "per_vertex_normals.h"
 #include "point_mesh_squared_distance.h"
+#include "pseudonormal_test.h"
 
 
 IGL_INLINE void igl::signed_distance(
@@ -117,6 +118,58 @@ IGL_INLINE double igl::signed_distance_pseudonormal(
   return s*sqrt(sqrd);
 }
 
+IGL_INLINE void igl::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)
+{
+  using namespace Eigen;
+  const size_t np = P.rows();
+  S.resize(np,1);
+  I.resize(np,1);
+  N.resize(np,3);
+  C.resize(np,3);
+# pragma omp parallel for if(np>1000)
+  for(size_t p = 0;p<np;p++)
+  {
+    double s,sqrd;
+    RowVector3d n,c;
+    int i = -1;
+    RowVector3d q = P.row(p);
+    signed_distance_pseudonormal(tree,V,F,FN,VN,EN,EMAP,q,s,sqrd,i,c,n);
+    S(p) = s*sqrt(sqrd);
+    I(p) = i;
+    N.row(p) = n;
+    C.row(p) = c;
+  }
+//  igl::AABB<MatrixXd,3> tree_P;
+//  MatrixXi J = VectorXi::LinSpaced(P.rows(),0,P.rows()-1);
+//  tree_P.init(P,J);
+//  tree.squared_distance(V,F,tree_P,P,J,S,I,C);
+//# pragma omp parallel for if(np>1000)
+//  for(size_t p = 0;p<np;p++)
+//  {
+//    RowVector3d c = C.row(p);
+//    RowVector3d q = P.row(p);
+//    const int f = I(p);
+//    double s;
+//    RowVector3d n;
+//    pseudonormal_test(V,F,FN,VN,EN,EMAP,q,f,c,s,n);
+//    N.row(p) = n;
+//    S(p) = s*sqrt(S(p));
+//  }
+
+}
+
 IGL_INLINE void igl::signed_distance_pseudonormal(
   const AABB<Eigen::MatrixXd,3> & tree,
   const Eigen::MatrixXd & V,
@@ -135,44 +188,7 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   using namespace Eigen;
   using namespace std;
   sqrd = tree.squared_distance(V,F,q,f,c);
-  const auto & qc = q-c;
-  RowVector3d b;
-  AABB<Eigen::MatrixXd,3>::barycentric_coordinates(
-    c,V.row(F(f,0)),V.row(F(f,1)),V.row(F(f,2)),b);
-  // Determine which normal to use
-  const double epsilon = 1e-12;
-  const int type = (b.array()<=epsilon).cast<int>().sum();
-  switch(type)
-  {
-    case 2:
-      // Find vertex
-      for(int x = 0;x<3;x++)
-      {
-        if(b(x)>epsilon)
-        {
-          n = VN.row(F(f,x));
-          break;
-        }
-      }
-      break;
-    case 1:
-      // Find edge
-      for(int x = 0;x<3;x++)
-      {
-        if(b(x)<=epsilon)
-        {
-          n = EN.row(EMAP(F.rows()*x+f));
-          break;
-        }
-      }
-      break;
-    default:
-      assert(false && "all barycentric coords zero.");
-    case 0:
-      n = FN.row(f);
-      break;
-  }
-  s = (qc.dot(n) >= 0 ? 1. : -1.);
+  pseudonormal_test(V,F,FN,VN,EN,EMAP,q,f,c,s,n);
 }
 
 IGL_INLINE double igl::signed_distance_winding_number(

include/igl/signed_distance.h

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_SIGNED_DISTANCE_H
 #define IGL_SIGNED_DISTANCE_H
+
 #include "igl_inline.h"
 #include "AABB.h"
 #include "WindingNumberAABB.h"
@@ -90,6 +91,19 @@ namespace igl
     int & i,
     Eigen::RowVector3d & c,
     Eigen::RowVector3d & n);
+  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);
 
   // Inputs:
   //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)

include/igl/triangle_triangle_adjacency.cpp

@@ -239,4 +239,5 @@ template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -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<int, -1, -1, 0, -1, -1> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -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<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(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<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/winding_number.cpp

@@ -213,4 +213,5 @@ IGL_INLINE void igl::winding_number_2(
 // Explicit template specialization
 template void igl::winding_number_2<double>(double const*, int, double const*, int, double const*, int, double*);
 template void igl::winding_number_3<double>(double const*, int, double const*, int, double const*, int, double*);
+template void igl::winding_number_3<double, int>(double const*, int, int const*, int, double const*, int, double*);
 #endif