general purpose AABB, replace cgal dependency for distance queries

Former-commit-id: c8a20e1ebde0805f5fe9b4bf7ae1d9289d444980

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.1.6  Major boolean robustness fix, drop CGAL dependency for AABB/distances
 1.1.5  Bug fix in booleans
 1.1.4  Edge collapsing and linear program solving
 1.1.3  Bug fixes in active set and boundary_conditions

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.1.5
+1.1.6

include/igl/AABB.h

@@ -65,13 +65,21 @@ namespace igl
       {
         swap(*this,other);
       }
-      ~AABB()
+      // Seems like there should have been an elegant solution to this using
+      // the copy-swap idiom above:
+      inline void deinit()
       {
+        m_element = -1;
+        m_box = Eigen::AlignedBox<Scalar,DIM>();
         delete m_left;
         m_left = NULL;
         delete m_right;
         m_right = NULL;
       }
+      ~AABB()
+      {
+        deinit();
+      }
       // Build an Axis-Aligned Bounding Box tree for a given mesh and given
       // serialization of a previous AABB tree.
       //
@@ -166,6 +174,7 @@ namespace igl
         const RowVectorDIMS & p,
         int & i,
         RowVectorDIMS & c) const;
+    private:
       inline Scalar squared_distance(
         const Eigen::PlainObjectBase<DerivedV> & V,
         const Eigen::MatrixXi & Ele, 
@@ -173,7 +182,7 @@ namespace igl
         const Scalar min_sqr_d,
         int & i,
         RowVectorDIMS & c) const;
-
+    public:
       template <
         typename DerivedP, 
         typename DerivedsqrD, 
@@ -203,6 +212,15 @@ namespace igl
         Scalar & sqr_d,
         int & i,
         RowVectorDIMS & c) const;
+    public:
+      template <int SS>
+        static
+        void barycentric_coordinates(
+          const RowVectorDIMS & p, 
+          const RowVectorDIMS & a, 
+          const RowVectorDIMS & b, 
+          const RowVectorDIMS & c,
+          Eigen::Matrix<Scalar,1,SS> & bary);
   };
 }
 
@@ -291,7 +309,7 @@ inline void igl::AABB<DerivedV,DIM>::init(
 {
   using namespace Eigen;
   using namespace std;
-  const int dim = V.cols();
+  assert(DIM == V.cols() && "V.cols() should matched declared dimension");
   const Scalar inf = numeric_limits<Scalar>::infinity();
   m_box = AlignedBox<Scalar,DIM>();
   // Compute bounding box
@@ -542,6 +560,8 @@ igl::AABB<DerivedV,DIM>::squared_distance(
   using namespace std;
   using namespace igl;
   Scalar sqr_d = min_sqr_d;
+  assert(DIM == 3 && "Code has only been tested for DIM == 3");
+  assert(Ele.cols() == 3 && "Code has only been tested for simplex size == 3");
 
   assert(m_element==-1 || (m_left == NULL && m_right == NULL));
   if(is_leaf())
@@ -646,30 +666,6 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
   using namespace igl;
   using namespace std;
 
-  // http://gamedev.stackexchange.com/a/23745
-  const auto & bary = [](
-      const RowVectorDIMS & p, 
-      const RowVectorDIMS & a, 
-      const RowVectorDIMS & b, 
-      const RowVectorDIMS & c, 
-      Scalar &u, 
-      Scalar &v, 
-      Scalar &w)
-  {
-    const RowVectorDIMS v0 = b - a;
-    const RowVectorDIMS v1 = c - a;
-    const RowVectorDIMS v2 = p - a;
-    Scalar d00 = v0.dot(v0);
-    Scalar d01 = v0.dot(v1);
-    Scalar d11 = v1.dot(v1);
-    Scalar d20 = v2.dot(v0);
-    Scalar d21 = v2.dot(v1);
-    Scalar denom = d00 * d11 - d01 * d01;
-    v = (d11 * d20 - d01 * d21) / denom;
-    w = (d00 * d21 - d01 * d20) / denom;
-    u = 1.0f - v - w;
-  };
-
   // Only one element per node
   // plane unit normal
   const RowVectorDIMS v10 = (V.row(Ele(m_element,1))- V.row(Ele(m_element,0)));
@@ -693,14 +689,14 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
     d_j = v.dot(un);
     pp = p - d_j*un;
     // determine if pp is inside triangle
-    Scalar b0,b1,b2;
-    bary(
-        pp,
-        V.row(Ele(m_element,0)),
-        V.row(Ele(m_element,1)),
-        V.row(Ele(m_element,2)),
-        b0,b1,b2);
-    inside_triangle = fabs(fabs(b0) + fabs(b1) + fabs(b2) - 1.) <= 1e-10;
+    Eigen::Matrix<Scalar,1,3> b;
+    barycentric_coordinates(
+          pp,
+          V.row(Ele(m_element,0)),
+          V.row(Ele(m_element,1)),
+          V.row(Ele(m_element,2)),
+          b);
+    inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
   }
   if(inside_triangle)
   {
@@ -754,4 +750,31 @@ inline void igl::AABB<DerivedV,DIM>::set_min(
   }
 }
 
+
+template <typename DerivedV, int DIM>
+template <int SS>
+inline void
+igl::AABB<DerivedV,DIM>::barycentric_coordinates(
+  const RowVectorDIMS & p, 
+  const RowVectorDIMS & a, 
+  const RowVectorDIMS & b, 
+  const RowVectorDIMS & c,
+  Eigen::Matrix<Scalar,1,SS> & bary)
+{
+  assert(SS==3);
+  // http://gamedev.stackexchange.com/a/23745
+  const RowVectorDIMS v0 = b - a;
+  const RowVectorDIMS v1 = c - a;
+  const RowVectorDIMS v2 = p - a;
+  Scalar d00 = v0.dot(v0);
+  Scalar d01 = v0.dot(v1);
+  Scalar d11 = v1.dot(v1);
+  Scalar d20 = v2.dot(v0);
+  Scalar d21 = v2.dot(v1);
+  Scalar denom = d00 * d11 - d01 * d01;
+  bary(1) = (d11 * d20 - d01 * d21) / denom;
+  bary(2) = (d00 * d21 - d01 * d20) / denom;
+  bary(0) = 1.0f - bary(1) - bary(2);
+};
+
 #endif

include/igl/InElementAABB.h

@@ -1,396 +0,0 @@
-#ifndef IGL_IN_ELEMENT_AABB_H
-#define IGL_IN_ELEMENT_AABB_H
-
-#include <Eigen/Core>
-#include <memory>
-#include <vector>
-namespace igl
-{
-  class InElementAABB
-  {
-    public:
-      std::shared_ptr<InElementAABB> m_left, m_right;
-      Eigen::RowVectorXd m_bb_min,m_bb_max;
-      // -1 non-leaf
-      int m_element;
-      InElementAABB():
-        m_left(NULL), m_right(NULL),
-        m_bb_min(), m_bb_max(), m_element(-1)
-      {}
-      // Build an Axis-Aligned Bounding Box tree for a given mesh and given
-      // serialization of a previous AABB tree.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. 
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
-      //   bb_mins  max_tree by dim list of bounding box min corner positions
-      //   bb_maxs  max_tree by dim list of bounding box max corner positions
-      //   elements  max_tree list of element or (not leaf id) indices into Ele
-      //   i  recursive call index {0}
-      inline void init(
-          const Eigen::MatrixXd & V,
-          const Eigen::MatrixXi & Ele, 
-          const Eigen::MatrixXd & bb_mins,
-          const Eigen::MatrixXd & bb_maxs,
-          const Eigen::VectorXi & elements,
-          const int i = 0);
-      // Wrapper for root with empty serialization
-      inline void init(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & Ele);
-      // Build an Axis-Aligned Bounding Box tree for a given mesh.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. 
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
-      //   SI  #Ele by dim list revealing for each coordinate where Ele's
-      //     barycenters would be sorted: SI(e,d) = i --> the dth coordinate of
-      //     the barycenter of the eth element would be placed at position i in a
-      //     sorted list.
-      //   I  #I list of indices into Ele of elements to include (for recursive
-      //     calls)
-      // 
-      inline void init(
-          const Eigen::MatrixXd & V,
-          const Eigen::MatrixXi & Ele, 
-          const Eigen::MatrixXi & SI,
-          const Eigen::VectorXi & I);
-      // Find the indices of elements containing given point.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. **Should be same as used to
-      //     construct mesh.**
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. **Should be same as used to
-      //     construct mesh.**
-      //   q  dim row-vector query position
-      //   first  whether to only return first element containing q
-      // Returns:
-      //   list of indices of elements containing q
-      inline std::vector<int> find(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & Ele, 
-        const Eigen::RowVectorXd & q, 
-        const bool first=false) const;
-  
-      // If number of elements m then total tree size should be 2*h where h is
-      // the deepest depth 2^ceil(log(#Ele*2-1))
-      inline int subtree_size();
-  
-      // Serialize this class into 3 arrays (so we can pass it pack to matlab)
-      //
-      // Outputs:
-      //   bb_mins  max_tree by dim list of bounding box min corner positions
-      //   bb_maxs  max_tree by dim list of bounding box max corner positions
-      //   elements  max_tree list of element or (not leaf id) indices into Ele
-      //   i  recursive call index into these arrays {0}
-      inline void serialize(
-        Eigen::MatrixXd & bb_mins,
-        Eigen::MatrixXd & bb_maxs,
-        Eigen::VectorXi & elements,
-        const int i = 0);
-  };
-}
-
-// Implementation
-#include <igl/volume.h>
-#include <igl/colon.h>
-#include <igl/doublearea.h>
-#include <igl/matlab_format.h>
-#include <igl/colon.h>
-#include <igl/sort.h>
-#include <igl/barycenter.h>
-#include <iostream>
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::MatrixXd & bb_mins,
-    const Eigen::MatrixXd & bb_maxs,
-    const Eigen::VectorXi & elements,
-    const int i)
-{
-  using namespace std;
-  using namespace Eigen;
-  if(bb_mins.size() > 0)
-  {
-    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
-    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
-    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
-    assert(bb_mins.rows() == elements.rows() &&
-      "Serial tree arrays must match");
-    // construct from serialization
-    m_bb_min = bb_mins.row(i);
-    m_bb_max = bb_maxs.row(i);
-    m_element = elements(i);
-    // Not leaf then recurse
-    if(m_element == -1)
-    {
-      m_left = make_shared<InElementAABB>();
-      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
-      m_right = make_shared<InElementAABB>();
-      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
-    }
-  }else
-  {
-    VectorXi allI = colon<int>(0,Ele.rows()-1);
-    MatrixXd BC;
-    barycenter(V,Ele,BC);
-    MatrixXi SI(BC.rows(),BC.cols());
-    {
-      MatrixXd _;
-      MatrixXi IS;
-      igl::sort(BC,1,true,_,IS);
-      // Need SI(i) to tell which place i would be sorted into
-      const int dim = IS.cols();
-      for(int i = 0;i<IS.rows();i++)
-      {
-        for(int d = 0;d<dim;d++)
-        {
-          SI(IS(i,d),d) = i;
-        }
-      }
-    }
-    init(V,Ele,SI,allI);
-  }
-}
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele)
-{
-  using namespace Eigen;
-  return init(V,Ele,MatrixXd(),MatrixXd(),VectorXi(),0);
-}
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::MatrixXi & SI,
-    const Eigen::VectorXi & I)
-{
-  using namespace Eigen;
-  using namespace std;
-  const int dim = V.cols();
-  const double inf = numeric_limits<double>::infinity();
-  m_bb_min.setConstant(1,dim, inf);
-  m_bb_max.setConstant(1,dim,-inf);
-  // Compute bounding box
-  for(int i = 0;i<I.rows();i++)
-  {
-    for(int c = 0;c<Ele.cols();c++)
-    {
-      for(int d = 0;d<dim;d++)
-      {
-        m_bb_min(d) = min(m_bb_min(d),V(Ele(I(i),c),d));
-        m_bb_max(d) = max(m_bb_max(d),V(Ele(I(i),c),d));
-      }
-    }
-  }
-  switch(I.size())
-  {
-    case 0:
-      {
-        assert(false);
-      }
-    case 1:
-      {
-        m_element = I(0);
-        break;
-      }
-    default:
-      {
-        // Compute longest direction
-        int max_d = -1;
-        double max_len = -inf;
-        for(int d = 0;d<dim;d++)
-        {
-          const auto diff = (m_bb_max[d] - m_bb_min[d]);
-          if( diff > max_len )
-          {
-            max_len = diff;
-            max_d = d;
-          }
-        }
-        // Can't use median on BC directly because many may have same value,
-        // but can use median on sorted BC indices
-        VectorXi SIdI(I.rows());
-        for(int i = 0;i<I.rows();i++)
-        {
-          SIdI(i) = SI(I(i),max_d);
-        }
-        // Since later I use <= I think I don't need to worry about odd/even
-        // Pass by copy to avoid changing input
-        const auto median = [](VectorXi A)->double
-        {
-          size_t n = A.size()/2;
-          nth_element(A.data(),A.data()+n,A.data()+A.size());
-          if(A.rows() % 2 == 1)
-          {
-            return A(n);
-          }else
-          {
-            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
-            return 0.5*(A(n)+A(n-1));
-          }
-        };
-        const double med = median(SIdI);
-        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
-        assert(LI.rows()+RI.rows() == I.rows());
-        // Distribute left and right
-        {
-          int li = 0;
-          int ri = 0;
-          for(int i = 0;i<I.rows();i++)
-          {
-            if(SIdI(i)<=med)
-            {
-              LI(li++) = I(i);
-            }else
-            {
-              RI(ri++) = I(i);
-            }
-          }
-        }
-        if(LI.rows()>0)
-        {
-          m_left = make_shared<InElementAABB>();
-          m_left->init(V,Ele,SI,LI);
-        }
-        if(RI.rows()>0)
-        {
-          m_right = make_shared<InElementAABB>();
-          m_right->init(V,Ele,SI,RI);
-        }
-      }
-  }
-}
-
-inline std::vector<int> igl::InElementAABB::find(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::RowVectorXd & q, 
-    const bool first) const
-{
-  using namespace std;
-  using namespace Eigen;
-  bool inside = true;
-  const int dim = m_bb_max.size();
-  assert(q.size() == m_bb_max.size());
-  const double epsilon = igl::EPS<double>();
-  // Check if outside bounding box
-  for(int d = 0;d<q.size()&&inside;d++)
-  {
-    inside &= (q(d)-m_bb_min(d))>=epsilon;
-    inside &= (m_bb_max(d)-q(d))>=epsilon;
-  }
-  if(!inside)
-  {
-    return std::vector<int>();
-  }
-  if(m_element != -1)
-  {
-    // Initialize to some value > -epsilon
-    double a1=1,a2=1,a3=1,a4=1;
-    switch(dim)
-    {
-      case 3:
-        {
-          // Barycentric coordinates
-          const RowVector3d V1 = V.row(Ele(m_element,0));
-          const RowVector3d V2 = V.row(Ele(m_element,1));
-          const RowVector3d V3 = V.row(Ele(m_element,2));
-          const RowVector3d V4 = V.row(Ele(m_element,3));
-          a1 = volume_single(V2,V4,V3,(RowVector3d)q);
-          a2 = volume_single(V1,V3,V4,(RowVector3d)q);
-          a3 = volume_single(V1,V4,V2,(RowVector3d)q);
-          a4 = volume_single(V1,V2,V3,(RowVector3d)q);
-          break;
-        }
-      case 2:
-        {
-          // Barycentric coordinates
-          const Vector2d V1 = V.row(Ele(m_element,0));
-          const Vector2d V2 = V.row(Ele(m_element,1));
-          const Vector2d V3 = V.row(Ele(m_element,2));
-          double a0 = doublearea_single(V1,V2,V3);
-          a1 = doublearea_single(V1,V2,(Vector2d)q);
-          a2 = doublearea_single(V2,V3,(Vector2d)q);
-          a3 = doublearea_single(V3,V1,(Vector2d)q);
-          break;
-        }
-      default:assert(false);
-    }
-    if(
-        a1>=-epsilon && 
-        a2>=-epsilon && 
-        a3>=-epsilon && 
-        a4>=-epsilon)
-    {
-      return std::vector<int>(1,m_element);
-    }else
-    {
-      return std::vector<int>();
-    }
-  }
-  std::vector<int> left = m_left->find(V,Ele,q,first);
-  if(first && !left.empty())
-  {
-    return left;
-  }
-  std::vector<int> right = m_right->find(V,Ele,q,first);
-  if(first)
-  {
-    return right;
-  }
-  left.insert(left.end(),right.begin(),right.end());
-  return left;
-}
-
-inline int igl::InElementAABB::subtree_size()
-{
-  // 1 for self
-  int n = 1;
-  int n_left = 0,n_right = 0;
-  if(m_left != NULL)
-  {
-    n_left = m_left->subtree_size();
-  }
-  if(m_right != NULL)
-  {
-    n_right = m_right->subtree_size();
-  }
-  n += 2*std::max(n_left,n_right);
-  return n;
-}
-
-
-inline void igl::InElementAABB::serialize(
-    Eigen::MatrixXd & bb_mins,
-    Eigen::MatrixXd & bb_maxs,
-    Eigen::VectorXi & elements,
-    const int i)
-{
-  using namespace std;
-  // Calling for root then resize output
-  if(i==0)
-  {
-    const int m = subtree_size();
-    //cout<<"m: "<<m<<endl;
-    bb_mins.resize(m,m_bb_min.size());
-    bb_maxs.resize(m,m_bb_max.size());
-    elements.resize(m,1);
-  }
-  //cout<<i<<" ";
-  bb_mins.row(i) = m_bb_min;
-  bb_maxs.row(i) = m_bb_max;
-  elements(i) = m_element;
-  if(m_left != NULL)
-  {
-    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
-  }
-  if(m_right != NULL)
-  {
-    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
-  }
-}
-#endif

include/igl/barycentric_coordinates.h

@@ -46,6 +46,8 @@ namespace igl
   // Outputs:
   //   L  #P by e list of barycentric coordinates
   //   
+  // Known bugs: this code is not tested (and probably will not work) for
+  // triangles and queries in 3D even if the query lives in/on the triangle.
   template <
     typename DerivedP,
     typename DerivedA,

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -9,7 +9,7 @@
 #include "mesh_to_cgal_triangle_list.h"
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance(
   const Eigen::MatrixXd & P,
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
@@ -30,7 +30,7 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 }
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance_precompute(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance_precompute(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   CGAL::AABB_tree<
@@ -78,7 +78,7 @@ IGL_INLINE void igl::point_mesh_squared_distance_precompute(
 }
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance(
   const Eigen::MatrixXd & P,
   const CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
@@ -120,8 +120,8 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
-template void igl::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cgal::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/cgal/point_mesh_squared_distance.h

@@ -5,75 +5,77 @@
 // 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_POINT_MESH_SQUARED_DISTANCE_H
-#define IGL_POINT_MESH_SQUARED_DISTANCE_H
+#ifndef IGL_CGAL_POINT_MESH_SQUARED_DISTANCE_H
+#define IGL_CGAL_POINT_MESH_SQUARED_DISTANCE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
 #include <vector>
 #include "CGAL_includes.hpp"
 namespace igl
 {
-  // Compute distances from a set of points P to a triangle mesh (V,F)
-  //
-  // 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
-  // Outputs:
-  //   sqrD  #P list of smallest squared distances
-  //   I  #P list of facet indices corresponding to smallest distances
-  //   C  #P by 3 list of closest points
-  //
-  // Known bugs: This only computes distances to triangles. So unreferenced
-  // vertices and degenerate triangles (segments) are ignored.
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance(
-    const Eigen::MatrixXd & P,
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    Eigen::VectorXd & sqrD,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C);
-  // Probably can do this in a way that we don't pass around `tree` and `T`
-  //
-  // Outputs:
-  //   tree  CGAL's AABB tree
-  //   T  list of CGAL triangles in order of F (for determining which was found
-  //     in computation)
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance_precompute(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+  namespace cgal
+  {
+    // Compute distances from a set of points P to a triangle mesh (V,F)
+    //
+    // 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
+    // Outputs:
+    //   sqrD  #P list of smallest squared distances
+    //   I  #P list of facet indices corresponding to smallest distances
+    //   C  #P by 3 list of closest points
+    //
+    // Known bugs: This only computes distances to triangles. So unreferenced
+    // vertices and degenerate triangles (segments) are ignored.
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance(
+      const Eigen::MatrixXd & P,
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      Eigen::VectorXd & sqrD,
+      Eigen::VectorXi & I,
+      Eigen::MatrixXd & C);
+    // Probably can do this in a way that we don't pass around `tree` and `T`
+    //
+    // Outputs:
+    //   tree  CGAL's AABB tree
+    //   T  list of CGAL triangles in order of F (for determining which was found
+    //     in computation)
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance_precompute(
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      CGAL::AABB_tree<
+        CGAL::AABB_traits<Kernel, 
+          CGAL::AABB_triangle_primitive<Kernel, 
+            typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+          >
         >
-      >
-    > & tree,
-    std::vector<CGAL::Triangle_3<Kernel> > & T);
-  // Inputs:
-  //  see above
-  // Outputs:
-  //  see above
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance(
-    const Eigen::MatrixXd & P,
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      > & tree,
+      std::vector<CGAL::Triangle_3<Kernel> > & T);
+    // Inputs:
+    //  see above
+    // Outputs:
+    //  see above
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance(
+      const Eigen::MatrixXd & P,
+      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,
-    Eigen::VectorXd & sqrD,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C);
-
+      > & tree,
+      const std::vector<CGAL::Triangle_3<Kernel> > & T,
+      Eigen::VectorXd & sqrD,
+      Eigen::VectorXi & I,
+      Eigen::MatrixXd & C);
+  }
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/signed_distance.cpp

@@ -1,292 +0,0 @@
-// 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 "signed_distance.h"
-#include "../per_vertex_normals.h"
-#include "../per_edge_normals.h"
-#include "../per_face_normals.h"
-#include "../get_seconds.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,FN,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<
-    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 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;
-  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;
-
-  pp = tree.closest_point_and_primitive(q);
-  Point_3 & p = pp.first;
-  const auto & qp = q-p;
-  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;
-  // barycentric coordinates
-  const auto & area = [&p,&t](const int i, const int j)->FT
-  {
-    return sqrt(Triangle_3(p,t.vertex(i),t.vertex(j)).squared_area());
-  };
-  Vector3d b(area(1,2),area(2,0),area(0,1));
-  b /= b.sum();
-  // 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 c = 0;c<3;c++)
-      {
-        if(b(c)>epsilon)
-        {
-          n = VN.row(F(f,c));
-          break;
-        }
-      }
-      break;
-    case 1:
-      // Find edge
-      for(int c = 0;c<3;c++)
-      {
-        if(b(c)<=epsilon)
-        {
-          n = EN.row(EMAP(F.rows()*c+f));
-          break;
-        }
-      }
-      break;
-    default:
-      assert(false && "all barycentric coords zero.");
-    case 0:
-      n = FN.row(f);
-      break;
-  }
-  s = (v.dot(n) >= 0 ? 1. : -1.);
-}
-
-template <typename Kernel>
-IGL_INLINE typename Kernel::FT 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 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;
-  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;
-  using namespace Eigen;
-  using namespace std;
-  using namespace igl;
-  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()); 
-  sqrd = qp.squared_length();
-  const double w = hier.winding_number(eq);
-  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
-template void igl::signed_distance<CGAL::Epick>( const Eigen::MatrixXd & , const Eigen::MatrixXd & , const Eigen::MatrixXi & , const SignedDistanceType , Eigen::VectorXd & , Eigen::VectorXi &, Eigen::MatrixXd & , Eigen::MatrixXd & );
-template CGAL::Epick::FT igl::signed_distance_winding_number<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Epick::Point_3 const&);
-template CGAL::Epick::FT igl::signed_distance_pseudonormal<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Epick::Point_3 const&);
-template void igl::signed_distance_pseudonormal<CGAL::Simple_cartesian<double> >(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Simple_cartesian<double>::Point_3 const&, CGAL::Simple_cartesian<double>::FT&, CGAL::Simple_cartesian<double>::FT&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > >::Point_and_primitive_id&, Eigen::Matrix<double, 3, 1, 0, 3, 1>&);
-template void igl::signed_distance_winding_number<CGAL::Simple_cartesian<double> >(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Simple_cartesian<double>::Point_3 const&, CGAL::Simple_cartesian<double>::FT&, CGAL::Simple_cartesian<double>::FT&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > >::Point_and_primitive_id&);
-#endif

include/igl/cgal/signed_distance_isosurface.cpp

@@ -8,21 +8,20 @@
 #include "signed_distance_isosurface.h"
 #include "point_mesh_squared_distance.h"
 #include "complex_to_mesh.h"
-#include "signed_distance.h"
 
-#include <igl/per_face_normals.h>
-#include <igl/per_edge_normals.h>
-#include <igl/per_vertex_normals.h>
-#include <igl/centroid.h>
-#include <igl/WindingNumberAABB.h>
-#include <igl/matlab_format.h>
-#include <igl/remove_unreferenced.h>
+#include "../AABB.h"
+#include "../per_face_normals.h"
+#include "../per_edge_normals.h"
+#include "../per_vertex_normals.h"
+#include "../centroid.h"
+#include "../WindingNumberAABB.h"
+#include "../matlab_format.h"
+#include "../remove_unreferenced.h"
 
 #include <CGAL/Surface_mesh_default_triangulation_3.h>
 #include <CGAL/Complex_2_in_triangulation_3.h>
 #include <CGAL/make_surface_mesh.h>
 #include <CGAL/Implicit_surface_3.h>
-#include <CGAL/Polyhedron_3.h>
 #include <CGAL/IO/output_surface_facets_to_polyhedron.h>
 // Axis-aligned bounding box tree for tet tri intersection
 #include <CGAL/AABB_tree.h>
@@ -42,6 +41,7 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   Eigen::MatrixXi & F)
 {
   using namespace std;
+  using namespace Eigen;
 
   // default triangulation for Surface_mesher
   typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
@@ -53,18 +53,9 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   typedef GT::FT FT;
   typedef std::function<FT (Point_3)> Function;
   typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
-  typedef CGAL::Polyhedron_3<GT> Polyhedron;
-  typedef GT::Kernel Kernel;
-  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-  typedef typename std::vector<Triangle_3>::iterator Iterator;
-  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Tree tree;
-  vector<Triangle_3 > T;
-  point_mesh_squared_distance_precompute(IV,IF,tree,T);
+  AABB<Eigen::MatrixXd,3> tree;
+  tree.init(IV,IF);
 
   Eigen::MatrixXd FN,VN,EN;
   Eigen::MatrixXi E;
@@ -108,36 +99,23 @@ IGL_INLINE bool igl::signed_distance_isosurface(
     case SIGNED_DISTANCE_TYPE_DEFAULT:
     case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
       fun = 
-        [&tree,&hier,&level](const Point_3 & q) -> FT
+        [&tree,&IV,&IF,&hier,&level](const Point_3 & q) -> FT
         {
-          return signed_distance_winding_number(tree,hier,q)-level;
+          const double sd = signed_distance_winding_number(
+            tree,IV,IF,hier,RowVector3d(q.x(),q.y(),q.z()));
+          return sd-level;
         };
       break;
     case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
-      fun = [&tree,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
+      fun = [&tree,&IV,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
         {
-          return 
-            igl::signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q) - 
-            level;
+          const double sd = 
+            igl::signed_distance_pseudonormal(
+              tree,IV,IF,FN,VN,EN,EMAP,RowVector3d(q.x(),q.y(),q.z()));
+          return sd- level;
         };
       break;
   }
-    //[&tree,&hier,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) ->FT
-    //{
-    //  const FT p = signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q);
-    //  const FT w = signed_distance_winding_number(tree,hier,q);
-    //  if(w*p < 0 && (fabs(w) > 0.1 || fabs(p) > 0.1))
-    //  {
-    //    cout<<"q=["<<q.x()<<","<<q.y()<<","<<q.z()<<"];"<<endl;
-    //    cout<<matlab_format(n.transpose().eval(),"n")<<endl;
-    //    cout<<matlab_format(b.transpose().eval(),"b")<<endl;
-    //    cout<<"Sig difference: "<<type<<endl;
-    //    cout<<"w: "<<w<<endl;
-    //    cout<<"p: "<<p<<endl;
-    //    exit(1);
-    //  }
-    //  return w;
-    //},
   Sphere_3 bounding_sphere(cmid, (r+level)*(r+level));
   Surface_3 surface(fun,bounding_sphere);
   CGAL::Surface_mesh_default_criteria_3<Tr> 

include/igl/cgal/signed_distance_isosurface.h

@@ -7,8 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_SIGNED_DISTANCE_ISOSURFACE_H
 #define IGL_SIGNED_DISTANCE_ISOSURFACE_H
-#include <igl/igl_inline.h>
-#include "signed_distance.h"
+#include "../igl_inline.h"
+#include "../signed_distance.h"
 #include <Eigen/Core>
 namespace igl
 {
@@ -22,7 +22,8 @@ namespace igl
   //   angle_bound  lower bound on triangle angles (mesh quality) (e.g. 28)
   //   radius_bound  upper bound on triangle size (mesh density?) (e.g. 0.02)
   //   distance_bound  cgal mysterious parameter (mesh density?) (e.g. 0.01)
-  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  //   sign_type  method for computing distance _sign_ (see
+  //     ../signed_distance.h)
   // Outputs:
   //   V  #V by 3 list of input mesh vertex positions
   //   F  #F by 3 list of input triangle indices

include/igl/in_element.cpp

@@ -1,17 +1,18 @@
 #include "in_element.h"
 
+template <typename DerivedV, typename DerivedQ, int DIM>
 IGL_INLINE void igl::in_element(
-  const Eigen::MatrixXd & V,
+  const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::MatrixXi & Ele,
-  const Eigen::MatrixXd & Q,
-  const InElementAABB & aabb,
+  const Eigen::PlainObjectBase<DerivedQ> & Q,
+  const AABB<DerivedV,DIM> & aabb,
   Eigen::VectorXi & I)
 {
   using namespace std;
   using namespace Eigen;
   const int Qr = Q.rows();
   I.setConstant(Qr,1,-1);
-#pragma omp parallel for
+#pragma omp parallel for if (Qr>10000)
   for(int e = 0;e<Qr;e++)
   {
     // find all
@@ -23,20 +24,21 @@ IGL_INLINE void igl::in_element(
   }
 }
 
+template <typename DerivedV, typename DerivedQ, int DIM, typename Scalar>
 IGL_INLINE void igl::in_element(
-  const Eigen::MatrixXd & V,
+  const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::MatrixXi & Ele,
-  const Eigen::MatrixXd & Q,
-  const InElementAABB & aabb,
-  Eigen::SparseMatrix<double> & I)
+  const Eigen::PlainObjectBase<DerivedQ> & Q,
+  const AABB<DerivedV,DIM> & aabb,
+  Eigen::SparseMatrix<Scalar> & I)
 {
   using namespace std;
   using namespace Eigen;
   using namespace igl;
   const int Qr = Q.rows();
-  std::vector<Triplet<double> > IJV;
+  std::vector<Triplet<Scalar> > IJV;
   IJV.reserve(Qr);
-#pragma omp parallel for
+#pragma omp parallel for if (Qr>10000)
   for(int e = 0;e<Qr;e++)
   {
     // find all
@@ -44,7 +46,7 @@ IGL_INLINE void igl::in_element(
     for(const auto r : R)
     {
 #pragma omp critical
-      IJV.push_back(Triplet<double>(e,r,1));
+      IJV.push_back(Triplet<Scalar>(e,r,1));
     }
   }
   I.resize(Qr,Ele.rows());

include/igl/in_element.h

@@ -1,8 +1,15 @@
+// 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_IN_ELEMENT_H
 #define IGL_IN_ELEMENT_H
 
-#include <igl/igl_inline.h>
-#include "InElementAABB.h"
+#include "igl_inline.h"
+#include "AABB.h"
 #include <Eigen/Core>
 
 namespace igl
@@ -10,6 +17,8 @@ namespace igl
   // Determine whether each point in a list of points is in the elements of a
   // mesh.
   //
+  // templates:
+  //   DIM  dimension of vertices in V (# of columns)
   // Inputs:
   //   V  #V by dim list of mesh vertex positions. 
   //   Ele  #Ele by dim+1 list of mesh indices into #V. 
@@ -18,20 +27,22 @@ namespace igl
   // Outputs:
   //   I  #Q list of indices into Ele of first containing element (-1 means no
   //     containing element)
+  template <typename DerivedV, typename DerivedQ, int DIM>
   IGL_INLINE void in_element(
-    const Eigen::MatrixXd & V,
+    const Eigen::PlainObjectBase<DerivedV> & V,
     const Eigen::MatrixXi & Ele,
-    const Eigen::MatrixXd & Q,
-    const InElementAABB & aabb,
+    const Eigen::PlainObjectBase<DerivedQ> & Q,
+    const AABB<DerivedV,DIM> & aabb,
     Eigen::VectorXi & I);
   // Outputs:
   //   I  #Q by #Ele sparse matrix revealing whether each element contains each
   //     point: I(q,e) means point q is in element e
+  template <typename DerivedV, typename DerivedQ, int DIM, typename Scalar>
   IGL_INLINE void in_element(
-    const Eigen::MatrixXd & V,
+    const Eigen::PlainObjectBase<DerivedV> & V,
     const Eigen::MatrixXi & Ele,
-    const Eigen::MatrixXd & Q,
-    const InElementAABB & aabb,
+    const Eigen::PlainObjectBase<DerivedQ> & Q,
+    const AABB<DerivedV,DIM> & aabb,
     Eigen::SparseMatrix<double> & I);
   //
   // Example:

include/igl/point_mesh_squared_distance.cpp

@@ -0,0 +1,32 @@
+// 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 "point_mesh_squared_distance.h"
+#include "AABB.h"
+
+template <
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::point_mesh_squared_distance(
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  using namespace std;
+  AABB<DerivedV,3> tree;
+  tree.init(V,Ele);
+  return tree.squared_distance(V,Ele,P,sqrD,I,C);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/point_mesh_squared_distance.h

@@ -0,0 +1,47 @@
+// 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/.
+#ifndef IGL_POINT_MESH_SQUARED_DISTANCE_H
+#define IGL_POINT_MESH_SQUARED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Compute distances from a set of points P to a triangle mesh (V,F)
+  //
+  // 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
+  // Outputs:
+  //   sqrD  #P list of smallest squared distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices and degenerate triangles (segments) are ignored.
+  template <
+    typename DerivedP,
+    typename DerivedV,
+    typename DerivedsqrD,
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void point_mesh_squared_distance(
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+    Eigen::PlainObjectBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedC> & C);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_mesh_squared_distance.cpp"
+#endif
+
+#endif

include/igl/signed_distance.cpp

@@ -0,0 +1,204 @@
+// 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 "signed_distance.h"
+#include "per_vertex_normals.h"
+#include "per_edge_normals.h"
+#include "per_face_normals.h"
+#include "get_seconds.h"
+#include "point_mesh_squared_distance.h"
+
+
+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");
+
+  // Prepare distance computation
+  AABB<MatrixXd,3> tree;
+  tree.init(V,F);
+
+  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,FN,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 RowVector3d q = P.row(p);
+    double s,sqrd;
+    RowVector3d c;
+    int i=-1;
+    switch(sign_type)
+    {
+      default:
+        assert(false && "Unknown SignedDistanceType");
+      case SIGNED_DISTANCE_TYPE_DEFAULT:
+      case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+        signed_distance_winding_number(tree,V,F,hier,q,s,sqrd,i,c);
+        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);
+        N.row(p) = n;
+        break;
+      }
+    }
+    I(p) = i;
+    S(p) = s*sqrt(sqrd);
+    C.row(p) = c;
+  }
+}
+
+
+IGL_INLINE double igl::signed_distance_pseudonormal(
+  const AABB<Eigen::MatrixXd,3> & tree,
+  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)
+{
+  double s,sqrd;
+  Eigen::RowVector3d n,c;
+  int i = -1;
+  signed_distance_pseudonormal(tree,V,F,FN,VN,EN,EMAP,q,s,sqrd,i,c,n);
+  return s*sqrt(sqrd);
+}
+
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const AABB<Eigen::MatrixXd,3> & tree,
+  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,
+  double & s,
+  double & sqrd,
+  int & f,
+  Eigen::RowVector3d & c,
+  Eigen::RowVector3d & n)
+{
+  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.);
+}
+
+IGL_INLINE double 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)
+{
+  double s,sqrd;
+  Eigen::RowVector3d c;
+  int i=-1;
+  signed_distance_winding_number(tree,V,F,hier,q,s,sqrd,i,c);
+  return s*sqrt(sqrd);
+}
+
+
+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,
+  double & s,
+  double & sqrd,
+  int & i,
+  Eigen::RowVector3d & c)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  sqrd = tree.squared_distance(V,F,q,i,c);
+  const double w = hier.winding_number(q.transpose());
+  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
+#endif

include/igl/cgal/signed_distance.h → include/igl/signed_distance.h

@@ -7,11 +7,11 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_SIGNED_DISTANCE_H
 #define IGL_SIGNED_DISTANCE_H
-#include <igl/igl_inline.h>
-#include <igl/WindingNumberAABB.h>
+#include "igl_inline.h"
+#include "AABB.h"
+#include "WindingNumberAABB.h"
 #include <Eigen/Core>
 #include <vector>
-#include "CGAL_includes.hpp"
 namespace igl
 {
   enum SignedDistanceType
@@ -24,9 +24,6 @@ namespace igl
   };
   // 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
@@ -41,7 +38,6 @@ namespace igl
   //
   // Known bugs: This only computes distances to triangles. So unreferenced
   // vertices and degenerate triangles are ignored.
-  template <typename Kernel>
   IGL_INLINE void signed_distance(
     const Eigen::MatrixXd & P,
     const Eigen::MatrixXd & V,
@@ -53,12 +49,8 @@ namespace igl
     Eigen::MatrixXd & N);
   // Computes signed distance to mesh
   //
-  // Templates:
-  //   Kernal  CGAL computation and construction kernel (e.g.
-  //     CGAL::Simple_cartesian<double>)
   // Inputs:
-  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
-  //   T  #F list of CGAL triangles (see point_mesh_squared_distance.h)
+  //   tree  AABB acceleration tree (see AABB.h)
   //   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)
@@ -67,92 +59,61 @@ namespace igl
   //   q  Query point
   // Returns signed distance to mesh
   //
-  template <typename Kernel>
-  IGL_INLINE typename Kernel::FT 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,
+  IGL_INLINE double signed_distance_pseudonormal(
+    const AABB<Eigen::MatrixXd,3> & tree,
+    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 typename Kernel::Point_3 & q);
+    const Eigen::RowVector3d & q);
   // Outputs:
   //   s  sign
   //   sqrd  squared distance
-  //   pp  closest point and primitve
+  //   i  closest primitive
+  //   c  closest point
   //   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 AABB<Eigen::MatrixXd,3> & tree,
+    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 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);
+    const Eigen::RowVector3d & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::RowVector3d & c,
+    Eigen::RowVector3d & n);
 
   // Inputs:
-  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
+  //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)
   //   hier  Winding number evaluation hierarchy
   //   q  Query point
   // Returns signed distance to mesh
-  template <typename Kernel>
-  IGL_INLINE typename Kernel::FT 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,
+  IGL_INLINE double 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 typename Kernel::Point_3 & q);
+    const Eigen::RowVector3d & 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 AABB<Eigen::MatrixXd,3> & tree,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
     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);
+    const Eigen::RowVector3d & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::RowVector3d & c);
 }
 
 #ifndef IGL_STATIC_LIBRARY