extract point-simplex distance and triangle barycentric coordinates from AABB, use MatrixBase for const input params

Former-commit-id: 954c36bf14aaa03fc92c4d39705fc3f6131c8acf

include/igl/AABB.h

@@ -284,14 +284,6 @@ private:
         Scalar & sqr_d,
         int & i,
         RowVectorDIMS & c) const;
-public:
-      static
-      inline void barycentric_coordinates(
-        const RowVectorDIMS & p, 
-        const RowVectorDIMS & a, 
-        const RowVectorDIMS & b, 
-        const RowVectorDIMS & c,
-        Eigen::Matrix<Scalar,1,3> & bary);
 public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
     };
@@ -300,10 +292,12 @@ public:
 // Implementation
 #include "EPS.h"
 #include "barycenter.h"
+#include "barycentric_coordinates.h"
 #include "colon.h"
 #include "colon.h"
 #include "doublearea.h"
 #include "matlab_format.h"
+#include "point_simplex_squared_distance.h"
 #include "project_to_line_segment.h"
 #include "sort.h"
 #include "volume.h"
@@ -1026,101 +1020,11 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
 {
   using namespace Eigen;
   using namespace std;
-
-  // Simplex size
-  const size_t ss = Ele.cols();
-  // Only one element per node
-  // plane unit normal
-  bool inside_triangle = false;
-  Scalar d_j = std::numeric_limits<Scalar>::infinity();
-  RowVectorDIMS pp;
-  // Only consider triangles, and non-degenerate triangles at that
-  if(ss == 3 && 
-      Ele(m_primitive,0) != Ele(m_primitive,1) && 
-      Ele(m_primitive,1) != Ele(m_primitive,2) && 
-      Ele(m_primitive,2) != Ele(m_primitive,0))
-  {
-    assert(DIM == 3 && "Only implemented for 3D triangles");
-    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
-    // can't be const because of annoying DIM template
-    RowVector3S v10(0,0,0);
-    v10.head(DIM) = (V.row(Ele(m_primitive,1))- V.row(Ele(m_primitive,0)));
-    RowVector3S v20(0,0,0);
-    v20.head(DIM) = (V.row(Ele(m_primitive,2))- V.row(Ele(m_primitive,0)));
-    const RowVectorDIMS n = (v10.cross(v20)).head(DIM);
-    Scalar n_norm = n.norm();
-    if(n_norm > 0)
-    {
-      const RowVectorDIMS un = n/n.norm();
-      // vector to plane
-      const RowVectorDIMS bc = 
-        1./3.*
-        ( V.row(Ele(m_primitive,0))+
-          V.row(Ele(m_primitive,1))+
-          V.row(Ele(m_primitive,2)));
-      const auto & v = p-bc;
-      // projected point on plane
-      d_j = v.dot(un);
-      pp = p - d_j*un;
-      // determine if pp is inside triangle
-      Eigen::Matrix<Scalar,1,3> b;
-      barycentric_coordinates(
-            pp,
-            V.row(Ele(m_primitive,0)),
-            V.row(Ele(m_primitive,1)),
-            V.row(Ele(m_primitive,2)),
-            b);
-      inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
-    }
-  }
-  const auto & point_point_squared_distance = [&](const RowVectorDIMS & s)
-  {
-    const Scalar sqr_d_s = (p-s).squaredNorm();
-    set_min(p,sqr_d_s,m_primitive,s,sqr_d,i,c);
-  };
-  if(inside_triangle)
-  {
-    // point-triangle squared distance
-    const Scalar sqr_d_j = d_j*d_j;
-    //cout<<"point-triangle..."<<endl;
-    set_min(p,sqr_d_j,m_primitive,pp,sqr_d,i,c);
-  }else
-  {
-    if(ss >= 2)
-    {
-      // point-segment distance
-      // number of edges
-      size_t ne = ss==3?3:1;
-      for(size_t x = 0;x<ne;x++)
-      {
-        const size_t e1 = Ele(m_primitive,(x+1)%ss);
-        const size_t e2 = Ele(m_primitive,(x+2)%ss);
-        const RowVectorDIMS & s = V.row(e1);
-        const RowVectorDIMS & d = V.row(e2);
-        // Degenerate edge
-        if(e1 == e2 || (s-d).squaredNorm()==0)
-        {
-          // only consider once
-          if(e1 < e2)
-          {
-            point_point_squared_distance(s);
-          }
-          continue;
-        }
-        Matrix<Scalar,1,1> sqr_d_j_x(1,1);
-        Matrix<Scalar,1,1> t(1,1);
-        project_to_line_segment(p,s,d,t,sqr_d_j_x);
-        const RowVectorDIMS q = s+t(0)*(d-s);
-        set_min(p,sqr_d_j_x(0),m_primitive,q,sqr_d,i,c);
-      }
-    }else
-    {
-      // So then Ele is just a list of points...
-      assert(ss == 1);
-      const RowVectorDIMS & s = V.row(Ele(m_primitive,0));
-      point_point_squared_distance(s);
-    }
-  }
+  RowVectorDIMS c_candidate;
+  Scalar sqr_d_candidate;
+  igl::point_simplex_squared_distance<DIM>(
+    p,V,Ele,m_primitive,sqr_d_candidate,c_candidate);
+  set_min(p,sqr_d_candidate,m_primitive,c_candidate,sqr_d,i,c);
 }
 
 
@@ -1153,30 +1057,6 @@ inline void igl::AABB<DerivedV,DIM>::set_min(
 }
 
 
-template <typename DerivedV, int DIM>
-inline void
-igl::AABB<DerivedV,DIM>::barycentric_coordinates(
-  const RowVectorDIMS & p, 
-  const RowVectorDIMS & a, 
-  const RowVectorDIMS & b, 
-  const RowVectorDIMS & c,
-  Eigen::Matrix<Scalar,1,3> & bary)
-{
-  // 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);
-}
-
 template <typename DerivedV, int DIM>
 inline bool 
 igl::AABB<DerivedV,DIM>::intersect_ray(
@@ -1201,7 +1081,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // Actually process elements
     assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
     // Cheesecake way of hitting element
-    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hits);
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
   }
   std::vector<igl::Hit> left_hits;
   std::vector<igl::Hit> right_hits;
@@ -1248,7 +1128,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
       assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
       igl::Hit leaf_hit;
       if(
-        ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive).eval(),leaf_hit)&&
+        ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive),leaf_hit)&&
         leaf_hit.t < hit.t)
       {
         hit = leaf_hit;
@@ -1302,7 +1182,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // Actually process elements
     assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
     // Cheesecake way of hitting element
-    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hit);
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hit);
   }
 
   // Doesn't seem like smartly choosing left before/after right makes a

include/igl/barycentric_coordinates.cpp

@@ -55,33 +55,55 @@ template <
   typename DerivedC,
   typename DerivedL>
 IGL_INLINE void igl::barycentric_coordinates(
-  const Eigen::PlainObjectBase<DerivedP> & P,
-  const Eigen::PlainObjectBase<DerivedA> & A,
-  const Eigen::PlainObjectBase<DerivedB> & B,
-  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::MatrixBase<DerivedP> & P,
+  const Eigen::MatrixBase<DerivedA> & A,
+  const Eigen::MatrixBase<DerivedB> & B,
+  const Eigen::MatrixBase<DerivedC> & C,
   Eigen::PlainObjectBase<DerivedL> & L)
 {
   using namespace Eigen;
-  assert(P.cols() == 2 && "query must be in 2d");
-  assert(A.cols() == 2 && "corners must be in 2d");
-  assert(B.cols() == 2 && "corners must be in 2d");
-  assert(C.cols() == 2 && "corners must be in 2d");
+#ifndef NDEBUG
+  const int DIM = P.cols();
+  assert(A.cols() == DIM && "corners must be in same dimension as query");
+  assert(B.cols() == DIM && "corners must be in same dimension as query");
+  assert(C.cols() == DIM && "corners must be in same dimension as query");
   assert(P.rows() == A.rows() && "Must have same number of queries as corners");
   assert(A.rows() == B.rows() && "Corners must be same size");
   assert(A.rows() == C.rows() && "Corners must be same size");
-  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
-    VectorXS;
-  // Total area
-  VectorXS dblA,LA,LB,LC;
-  doublearea(P,B,C,LA);
-  doublearea(A,P,C,LB);
-  doublearea(A,B,P,LC);
-  doublearea(A,B,C,dblA);
+#endif
+
+  // http://gamedev.stackexchange.com/a/23745
+  typedef 
+    Eigen::Array<
+      typename DerivedP::Scalar,
+               DerivedP::RowsAtCompileTime,
+               DerivedP::ColsAtCompileTime>
+    ArrayS;
+  typedef 
+    Eigen::Array<
+      typename DerivedP::Scalar,
+               DerivedP::RowsAtCompileTime,
+               1>
+    VectorS;
+
+  const ArrayS v0 = B.array() - A.array();
+  const ArrayS v1 = C.array() - A.array();
+  const ArrayS v2 = P.array() - A.array();
+  VectorS d00 = (v0*v0).rowwise().sum();
+  VectorS d01 = (v0*v1).rowwise().sum();
+  VectorS d11 = (v1*v1).rowwise().sum();
+  VectorS d20 = (v2*v0).rowwise().sum();
+  VectorS d21 = (v2*v1).rowwise().sum();
+  VectorS denom = d00 * d11 - d01 * d01;
   L.resize(P.rows(),3);
-  L<<LA,LB,LC;
-  L.array().colwise() /= dblA.array();
+  L.col(1) = (d11 * d20 - d01 * d21) / denom;
+  L.col(2) = (d00 * d21 - d01 * d20) / denom;
+  L.col(0) = 1.0f -(L.col(1) + L.col(2)).array();
 }
 
 #ifdef IGL_STATIC_LIBRARY
 template void igl::barycentric_coordinates<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>, Eigen::Matrix<double, -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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 #endif

include/igl/barycentric_coordinates.h

@@ -39,15 +39,13 @@ namespace igl
   // Compute barycentric coordinates in a triangle
   //
   // Inputs:
-  //   P  #P by 2 Query points in 2d
-  //   A  #P by 2 Triangle corners in 2d
-  //   B  #P by 2 Triangle corners in 2d
-  //   C  #P by 2 Triangle corners in 2d
+  //   P  #P by dim Query points
+  //   A  #P by dim Triangle corners
+  //   B  #P by dim Triangle corners
+  //   C  #P by dim Triangle corners
   // Outputs:
-  //   L  #P by e list of barycentric coordinates
+  //   L  #P by 3 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,
@@ -55,10 +53,10 @@ namespace igl
     typename DerivedC,
     typename DerivedL>
   IGL_INLINE void barycentric_coordinates(
-    const Eigen::PlainObjectBase<DerivedP> & P,
-    const Eigen::PlainObjectBase<DerivedA> & A,
-    const Eigen::PlainObjectBase<DerivedB> & B,
-    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const Eigen::MatrixBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedL> & L);
 
 }

include/igl/copyleft/cgal/piecewise_constant_winding_number.cpp

@@ -24,6 +24,6 @@ IGL_INLINE bool igl::copyleft::cgal::piecewise_constant_winding_number(
   // resolve intersections
   remesh_self_intersections(V,F,{},VV,FF,IF,J,IM);
   // _apply_ duplicate vertex mapping IM to FF
-  for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
-  return piecewise_constant_winding_number(FF);
+  std::for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
+  return igl::piecewise_constant_winding_number(FF);
 }

include/igl/point_simplex_squared_distance.cpp

@@ -0,0 +1,148 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_simplex_squared_distance.h"
+#include "project_to_line_segment.h"
+#include "barycentric_coordinates.h"
+#include <Eigen/Geometry>
+#include <limits>
+#include <cassert>
+
+template <
+  int DIM,
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedsqr_d,
+  typename Derivedc>
+IGL_INLINE void igl::point_simplex_squared_distance(
+  const Eigen::PlainObjectBase<Derivedp> & p,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedEle> & Ele,
+  const typename DerivedEle::Index primitive,
+  Derivedsqr_d & sqr_d,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+  assert(p.size() == DIM);
+  assert(V.cols() == DIM);
+  assert(Ele.cols() <= DIM+1);
+  assert(Ele.cols() <= 3 && "Only simplices up to triangles are considered");
+
+  typedef Derivedsqr_d Scalar;
+  typedef typename Eigen::Matrix<Scalar,1,DIM> RowVectorDIMS;
+  sqr_d = std::numeric_limits<Scalar>::infinity();
+  const auto & set_min = 
+    [&sqr_d,&c](const Scalar & sqr_d_candidate, const RowVectorDIMS & c_candidate)
+  {
+    if(sqr_d_candidate < sqr_d)
+    {
+      sqr_d = sqr_d_candidate;
+      c = c_candidate;
+    }
+  };
+
+  // Simplex size
+  const size_t ss = Ele.cols();
+  // Only one element per node
+  // plane unit normal
+  bool inside_triangle = false;
+  Scalar d_j = std::numeric_limits<Scalar>::infinity();
+  RowVectorDIMS pp;
+  // Only consider triangles, and non-degenerate triangles at that
+  if(ss == 3 && 
+      Ele(primitive,0) != Ele(primitive,1) && 
+      Ele(primitive,1) != Ele(primitive,2) && 
+      Ele(primitive,2) != Ele(primitive,0))
+  {
+    assert(DIM == 3 && "Only implemented for 3D triangles");
+    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+    // can't be const because of annoying DIM template
+    RowVector3S v10(0,0,0);
+    v10.head(DIM) = (V.row(Ele(primitive,1))- V.row(Ele(primitive,0)));
+    RowVector3S v20(0,0,0);
+    v20.head(DIM) = (V.row(Ele(primitive,2))- V.row(Ele(primitive,0)));
+    const RowVectorDIMS n = (v10.cross(v20)).head(DIM);
+    Scalar n_norm = n.norm();
+    if(n_norm > 0)
+    {
+      const RowVectorDIMS un = n/n.norm();
+      // vector to plane
+      const RowVectorDIMS bc = 
+        1./3.*
+        ( V.row(Ele(primitive,0))+
+          V.row(Ele(primitive,1))+
+          V.row(Ele(primitive,2)));
+      const auto & v = p-bc;
+      // projected point on plane
+      d_j = v.dot(un);
+      pp = p - d_j*un;
+      // determine if pp is inside triangle
+      Eigen::Matrix<Scalar,1,3> b;
+      barycentric_coordinates(
+            pp,
+            V.row(Ele(primitive,0)),
+            V.row(Ele(primitive,1)),
+            V.row(Ele(primitive,2)),
+            b);
+      inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
+    }
+  }
+  const auto & point_point_squared_distance = [&](const RowVectorDIMS & s)
+  {
+    const Scalar sqr_d_s = (p-s).squaredNorm();
+    set_min(sqr_d_s,s);
+  };
+  if(inside_triangle)
+  {
+    // point-triangle squared distance
+    const Scalar sqr_d_j = d_j*d_j;
+    //cout<<"point-triangle..."<<endl;
+    set_min(sqr_d_j,pp);
+  }else
+  {
+    if(ss >= 2)
+    {
+      // point-segment distance
+      // number of edges
+      size_t ne = ss==3?3:1;
+      for(size_t x = 0;x<ne;x++)
+      {
+        const size_t e1 = Ele(primitive,(x+1)%ss);
+        const size_t e2 = Ele(primitive,(x+2)%ss);
+        const RowVectorDIMS & s = V.row(e1);
+        const RowVectorDIMS & d = V.row(e2);
+        // Degenerate edge
+        if(e1 == e2 || (s-d).squaredNorm()==0)
+        {
+          // only consider once
+          if(e1 < e2)
+          {
+            point_point_squared_distance(s);
+          }
+          continue;
+        }
+        Eigen::Matrix<Scalar,1,1> sqr_d_j_x(1,1);
+        Eigen::Matrix<Scalar,1,1> t(1,1);
+        project_to_line_segment(p,s,d,t,sqr_d_j_x);
+        const RowVectorDIMS q = s+t(0)*(d-s);
+        set_min(sqr_d_j_x(0),q);
+      }
+    }else
+    {
+      // So then Ele is just a list of points...
+      assert(ss == 1);
+      const RowVectorDIMS & s = V.row(Ele(primitive,0));
+      point_point_squared_distance(s);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+#endif

include/igl/point_simplex_squared_distance.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_SIMPLEX_SQUARED_DISTANCE_H
+#define IGL_POINT_SIMPLEX_SQUARED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Determine squared distance from a point to linear simplex
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::PlainObjectBase<Derivedp> & p,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_simplex_squared_distance.cpp"
+#endif
+#endif

include/igl/pseudonormal_test.cpp

@@ -6,7 +6,8 @@
 // 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 "barycentric_coordinates.h"
+#include "doublearea.h"
 #include "project_to_line_segment.h"
 #include <cassert>
 
@@ -49,7 +50,7 @@ IGL_INLINE void igl::pseudonormal_test(
   const double epsilon = 1e-12;
   if(area>MIN_DOUBLE_AREA)
   {
-    AABB<Eigen::MatrixXd,3>::barycentric_coordinates( c,A,B,C,b);
+    barycentric_coordinates( c,A,B,C,b);
     // Determine which normal to use
     const int type = (b.array()<=epsilon).cast<int>().sum();
     switch(type)

include/igl/ray_mesh_intersect.cpp

@@ -11,10 +11,10 @@ template <
   typename DerivedV, 
   typename DerivedF> 
 IGL_INLINE bool igl::ray_mesh_intersect(
-  const Eigen::PlainObjectBase<Derivedsource> & s,
-  const Eigen::PlainObjectBase<Deriveddir> & dir,
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedsource> & s,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   std::vector<igl::Hit> & hits)
 {
   using namespace Eigen;
@@ -53,10 +53,10 @@ template <
   typename DerivedV, 
   typename DerivedF> 
 IGL_INLINE bool igl::ray_mesh_intersect(
-  const Eigen::PlainObjectBase<Derivedsource> & source,
-  const Eigen::PlainObjectBase<Deriveddir> & dir,
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedsource> & source,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   igl::Hit & hit)
 {
   std::vector<igl::Hit> hits;

include/igl/ray_mesh_intersect.h

@@ -23,10 +23,10 @@ namespace igl
     typename DerivedV, 
     typename DerivedF> 
   IGL_INLINE bool ray_mesh_intersect(
-    const Eigen::PlainObjectBase<Derivedsource> & source,
-    const Eigen::PlainObjectBase<Deriveddir> & dir,
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & dir,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
     std::vector<igl::Hit> & hits);
   // Outputs:
   //   hit  first hit, set only if it exists
@@ -37,10 +37,10 @@ namespace igl
     typename DerivedV, 
     typename DerivedF> 
   IGL_INLINE bool ray_mesh_intersect(
-    const Eigen::PlainObjectBase<Derivedsource> & source,
-    const Eigen::PlainObjectBase<Deriveddir> & dir,
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & dir,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
     igl::Hit & hit);
 }
 #ifndef IGL_STATIC_LIBRARY