Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 61f777753f0d2c2f1a7b25f1cbfaf3f9b1541099

include/igl/AABB.h

@@ -768,7 +768,7 @@ inline typename igl::AABB<DerivedV,DIM>::Scalar igl::AABB<DerivedV,DIM>::squared
   const AABB<Derivedother_V,DIM> * other,
   const Eigen::PlainObjectBase<Derivedother_V> & other_V,
   const Eigen::MatrixXi & other_Ele, 
-  const Scalar min_sqr_d,
+  const Scalar /*min_sqr_d*/,
   Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
   Eigen::PlainObjectBase<DerivedI> & I,
   Eigen::PlainObjectBase<DerivedC> & C) const
@@ -1098,7 +1098,11 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
 
 template <typename DerivedV, int DIM>
 inline void igl::AABB<DerivedV,DIM>::set_min(
-  const RowVectorDIMS & p,
+  const RowVectorDIMS & 
+#ifndef NDEBUG
+  p
+#endif
+  ,
   const Scalar sqr_d_candidate,
   const int i_candidate,
   const RowVectorDIMS & c_candidate,

include/igl/WindingNumberTree.h

@@ -406,13 +406,16 @@ inline void igl::WindingNumberTree<Point>::print(const char * tab)
 }
 
 template <typename Point>
-inline double igl::WindingNumberTree<Point>::max_abs_winding_number(const Point & p) const
+inline double 
+igl::WindingNumberTree<Point>::max_abs_winding_number(const Point & /*p*/) const
 {
   return std::numeric_limits<double>::infinity();
 }
 
 template <typename Point>
-inline double igl::WindingNumberTree<Point>::max_simple_abs_winding_number(const Point & p) const
+inline double 
+igl::WindingNumberTree<Point>::max_simple_abs_winding_number(
+  const Point & /*p*/) const
 {
   using namespace std;
   return numeric_limits<double>::infinity();

include/igl/angle_bound_frame_fields.h

@@ -6,8 +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/.
 
-#ifndef IGL_ANGLE_BOUND_FRAME_FIELDS
-#define IGL_ANGLE_BOUND_FRAME_FIELDS
+#ifndef IGL_ANGLE_BOUND_FRAME_FIELDS_H
+#define IGL_ANGLE_BOUND_FRAME_FIELDS_H
 #include "igl_inline.h"
 
 #include <Eigen/Core>
@@ -58,4 +58,4 @@ namespace igl {
 #endif
 
 
-#endif /* defined(IGL_ANGLE_BOUND_FRAME_FIELDS) */
+#endif 

include/igl/arap_dof.cpp

@@ -184,7 +184,7 @@ IGL_INLINE bool igl::arap_dof_precomputation(
     //printf("CSM_M(): Mi\n");
     LbsMatrixType M_i;
     //printf("CSM_M(): slice\n");
-    slice(M,(span_n.array()+i*n).matrix(),span_mlbs_cols,M_i);
+    slice(M,(span_n.array()+i*n).matrix().eval(),span_mlbs_cols,M_i);
     LbsMatrixType M_i_dim;
     data.CSM_M[i].resize(k*data.dim,data.m*data.dim*(data.dim+1));
     assert(data.CSM_M[i].cols() == M.cols());

include/igl/boolean/CSGTree.h

@@ -92,7 +92,10 @@ namespace igl
           const MeshBooleanType & type)
         {
           // conduct boolean operation
-          mesh_boolean(A.V(),A.F(),B.V(),B.F(),type,m_V,m_F,m_J);
+          {
+            Eigen::VectorXi I;
+            mesh_boolean(A.V(),A.F(),B.V(),B.F(),type,m_V,m_F,m_J,I);
+          }
           // reindex m_J
           std::for_each(m_J.data(),m_J.data()+m_J.size(),
             [&](typename VectorJ::Scalar & j)

include/igl/boolean/mesh_boolean.cpp

@@ -19,62 +19,6 @@
 #include <iostream>
 
 //#define IGL_MESH_BOOLEAN_DEBUG
-
-template <
-  typename DerivedVA,
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename DerivedVC,
-  typename DerivedFC,
-  typename DerivedJ>
-IGL_INLINE void igl::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC,
-  Eigen::PlainObjectBase<DerivedJ > & J)
-{
-  const std::function<void(
-    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)>
-    empty_fun;
-  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
-}
-
-template <
-  typename DerivedVA,
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename DerivedVC,
-  typename DerivedFC>
-IGL_INLINE void igl::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC)
-{
-  Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
-  const std::function<void(
-    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&)>
-    empty_fun;
-  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
-}
-
 template <
   typename DerivedVA,
   typename DerivedFA,
@@ -82,7 +26,8 @@ template <
   typename DerivedFB,
   typename DerivedVC,
   typename DerivedFC,
-  typename DerivedJ>
+  typename DerivedJ,
+  typename DerivedI>
 IGL_INLINE void igl::boolean::mesh_boolean(
   const Eigen::PlainObjectBase<DerivedVA > & VA,
   const Eigen::PlainObjectBase<DerivedFA > & FA,
@@ -92,13 +37,15 @@ IGL_INLINE void igl::boolean::mesh_boolean(
   const std::function<void(
     const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
     const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)>
-    & resolve_fun,
+    Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
+    Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+    Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&,
+    Eigen::Matrix<typename DerivedI::Scalar, Eigen::Dynamic,1>&)> &
+    resolve_fun,
   Eigen::PlainObjectBase<DerivedVC > & VC,
   Eigen::PlainObjectBase<DerivedFC > & FC,
-  Eigen::PlainObjectBase<DerivedJ > & J)
+  Eigen::PlainObjectBase<DerivedJ > & J,
+  Eigen::PlainObjectBase<DerivedI > & I)
 {
   using namespace Eigen;
   using namespace std;
@@ -114,7 +61,7 @@ IGL_INLINE void igl::boolean::mesh_boolean(
   typedef Matrix<ExactScalar,Dynamic,3> MatrixX3ES;
   typedef Matrix<Index,Dynamic,3> MatrixX3I;
   typedef Matrix<Index,Dynamic,2> MatrixX2I;
-  typedef Matrix<Index,Dynamic,1> VectorXI;
+  typedef Matrix<typename DerivedI::Scalar,Dynamic,1> VectorXI;
   typedef Matrix<typename DerivedJ::Scalar,Dynamic,1> VectorXJ;
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"mesh boolean..."<<endl;
@@ -149,18 +96,26 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     const MatrixX3I & F,
     MatrixX3ES & CV,
     MatrixX3I & CF,
-    VectorXJ & J)
+    VectorXJ & J,
+    VectorXI & I)
   {
     MatrixX3ES SV;
     MatrixX3I SF;
     MatrixX2I SIF;
-    VectorXI SIM,UIM;
     igl::cgal::RemeshSelfIntersectionsParam params;
-    remesh_self_intersections(V,F,params,SV,SF,SIF,J,SIM);
-    for_each(SF.data(),SF.data()+SF.size(),[&SIM](int & a){a=SIM(a);});
+    remesh_self_intersections(V,F,params,SV,SF,SIF,J,I);
+    for_each(SF.data(),SF.data()+SF.size(),
+      [&I](typename MatrixX3I::Scalar & a){a=I(a);});
     {
+      Eigen::Matrix<typename MatrixX3S::Index,Dynamic,1> UIM;
       remove_unreferenced(SV,SF,CV,CF,UIM);
+      for_each(I.data(),I.data()+I.size(),
+        [&UIM](typename VectorXI::Scalar & a){a=UIM(a);});
     }
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+    cout<<"#F:  "<<F.rows()<<endl;
+    cout<<"#CF: "<<CF.rows()<<endl;
+#endif
   };
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
@@ -170,12 +125,13 @@ IGL_INLINE void igl::boolean::mesh_boolean(
   MatrixX3ES EV;
   MatrixX3I CF;
   VectorXJ CJ;
+  Eigen::Matrix<typename DerivedI::Scalar,Dynamic, 1> CI;
   if(resolve_fun)
   {
-    resolve_fun(V,F,CV,CF,CJ);
+    resolve_fun(V,F,CV,CF,CJ,CI);
   }else
   {
-    libigl_resolve(V,F,EV,CF,CJ);
+    libigl_resolve(V,F,EV,CF,CJ,CI);
     CV.resize(EV.rows(), EV.cols());
     // Just use f'ing for loops. What if EV and CV don't use the same ordering?
     for(int i=0;i<EV.rows();i++)
@@ -192,6 +148,7 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     FC = CF;
     VC = CV;
     J = CJ;
+    I = CI;
     return;
   }
 
@@ -200,13 +157,13 @@ IGL_INLINE void igl::boolean::mesh_boolean(
 #endif
   Matrix<bool,Dynamic,1> from_A(CF.rows());
   // peel layers keeping track of odd and even flips
-  VectorXi I;
+  VectorXi iter;
   Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(EV,CF,I,flip);
+  peel_outer_hull_layers(EV,CF,iter,flip);
   //Array<bool,Dynamic,1> even = igl::mod(I,2).array()==0;
   const auto even = [&](const Index & f)->bool
   {
-    return (I(f)%2)==0;
+    return (iter(f)%2)==0;
   };
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
@@ -335,8 +292,15 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     }
   }
   // remove unreferenced vertices
-  VectorXi newIM;
-  remove_unreferenced(CV,G,VC,FC,newIM);
+  {
+    VectorXi newIM;
+    remove_unreferenced(CV,G,VC,FC,newIM);
+    I.resize(CI.rows(),CI.cols());
+    for(int i = 0;i<CI.rows();i++)
+    {
+      I(i) = newIM(CI(i));
+    }
+  }
   //cerr<<"warning not removing unref"<<endl;
   //VC = CV;
   //FC = G;
@@ -350,7 +314,73 @@ IGL_INLINE void igl::boolean::mesh_boolean(
 #endif
 }
 
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC,
+  typename DerivedJ,
+  typename DerivedI>
+IGL_INLINE void igl::boolean::mesh_boolean(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC,
+  Eigen::PlainObjectBase<DerivedJ > & J,
+  Eigen::PlainObjectBase<DerivedI > & I)
+{
+  const std::function<void(
+    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&,
+          Eigen::Matrix<typename DerivedI::Scalar, Eigen::Dynamic,1>&) >
+    empty_fun;
+  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J,I);
+}
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC>
+IGL_INLINE void igl::boolean::mesh_boolean(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC)
+{
+  Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
+  typedef Eigen::Matrix<typename DerivedVC::Index, Eigen::Dynamic,1> VectorXI;
+  VectorXI I;
+  const std::function<void(
+    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&,
+          Eigen::Matrix<typename VectorXI::Scalar, Eigen::Dynamic,1>&)>
+    empty_fun;
+  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J,I);
+}
+
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType 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<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::boolean::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // This is a hack to discuss. I'm not sure why this _doesn't_ create
 // duplicate symbols.
 #include <igl/remove_unreferenced.cpp>
@@ -361,16 +391,14 @@ igl::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>,
 #include <igl/cgal/outer_hull.cpp>
 template void igl::cgal::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 #include <igl/slice.cpp>
-template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
 #include <igl/barycenter.cpp>
 template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
 #include <igl/mod.cpp>
 template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > igl::mod<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int);
 #include <igl/outer_element.cpp>
 template void igl::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-// Explicit template specialization
-template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType 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<int, -1, 1, 0, -1, 1> >&);
-template void igl::boolean::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+#include <igl/colon.cpp>
+template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
+
 
 #endif

include/igl/boolean/mesh_boolean.h

@@ -21,82 +21,85 @@ namespace igl
     //  oriented meshes.
     // 
     //  Inputs:
-    //    V  #V by 3 list of vertex positions of first mesh
-    //    F  #F by 3 list of triangle indices into V
-    //    U  #U by 3 list of vertex positions of second mesh
-    //    G  #G by 3 list of triangle indices into U
+    //    VA  #VA by 3 list of vertex positions of first mesh
+    //    FA  #FA by 3 list of triangle indices into VA
+    //    VB  #VB by 3 list of vertex positions of second mesh
+    //    FB  #FB by 3 list of triangle indices into VB
     //    type  type of boolean operation
+    //    resolve_fun  function handle for computing resolve of a
+    //      self-intersections of a mesh and outputting the new mesh.
     //  Outputs:
-    //    W  #W by 3 list of vertex positions of boolean result mesh
-    //    H  #H by 3 list of triangle indices into W
-    //    J  #H list of indices into [FA;FB] revealing "birth" facet
-    //  
-    //  See also: self_intersect
+    //    VC  #VC by 3 list of vertex positions of boolean result mesh
+    //    FC  #FC by 3 list of triangle indices into VC
+    //    J  #FC list of indices into [FA;FB] revealing "birth" facet
+    //    I  #VA+#VB list of indices into SV (SV=[VA;VB;SVA;SVB], where SVA and
+    //      SVB are the new vertices from resolving intersections) revealing
+    //      "birth" vertices.
+    //
+    //  See also: mesh_boolean_cork, intersect_other, remesh_self_intersections
     //     
     template <
       typename DerivedVA,
-               typename DerivedFA,
-               typename DerivedVB,
-               typename DerivedFB,
-               typename DerivedVC,
-               typename DerivedFC,
-               typename DerivedJ>
-                 IGL_INLINE void mesh_boolean(
-                     const Eigen::PlainObjectBase<DerivedVA > & VA,
-                     const Eigen::PlainObjectBase<DerivedFA > & FA,
-                     const Eigen::PlainObjectBase<DerivedVB > & VB,
-                     const Eigen::PlainObjectBase<DerivedFB > & FB,
-                     const MeshBooleanType & type,
-                     Eigen::PlainObjectBase<DerivedVC > & VC,
-                     Eigen::PlainObjectBase<DerivedFC > & FC,
-                     Eigen::PlainObjectBase<DerivedJ > & J);
+      typename DerivedFA,
+      typename DerivedVB,
+      typename DerivedFB,
+      typename DerivedVC,
+      typename DerivedFC,
+      typename DerivedJ,
+      typename DerivedI>
+    IGL_INLINE void mesh_boolean(
+      const Eigen::PlainObjectBase<DerivedVA > & VA,
+      const Eigen::PlainObjectBase<DerivedFA > & FA,
+      const Eigen::PlainObjectBase<DerivedVB > & VB,
+      const Eigen::PlainObjectBase<DerivedFB > & FB,
+      const MeshBooleanType & type,
+      const std::function<void(
+        const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+        const Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
+        Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+        Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
+        Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1>&,
+        Eigen::Matrix<typename DerivedI::Scalar,Eigen::Dynamic,1>&)> &
+        resolve_fun,
+      Eigen::PlainObjectBase<DerivedVC > & VC,
+      Eigen::PlainObjectBase<DerivedFC > & FC,
+      Eigen::PlainObjectBase<DerivedJ > & J,
+      Eigen::PlainObjectBase<DerivedI > & I);
     template <
       typename DerivedVA,
-               typename DerivedFA,
-               typename DerivedVB,
-               typename DerivedFB,
-               typename DerivedVC,
-               typename DerivedFC>
-                 IGL_INLINE void mesh_boolean(
-                     const Eigen::PlainObjectBase<DerivedVA > & VA,
-                     const Eigen::PlainObjectBase<DerivedFA > & FA,
-                     const Eigen::PlainObjectBase<DerivedVB > & VB,
-                     const Eigen::PlainObjectBase<DerivedFB > & FB,
-                     const MeshBooleanType & type,
-                     Eigen::PlainObjectBase<DerivedVC > & VC,
-                     Eigen::PlainObjectBase<DerivedFC > & FC);
-    // Inputs:
-    //   resolve_fun  function handle for computing resolve of a
-    //     self-intersections of a mesh and outputting the new mesh.
+      typename DerivedFA,
+      typename DerivedVB,
+      typename DerivedFB,
+      typename DerivedVC,
+      typename DerivedFC,
+      typename DerivedJ,
+      typename DerivedI>
+    IGL_INLINE void mesh_boolean(
+      const Eigen::PlainObjectBase<DerivedVA > & VA,
+      const Eigen::PlainObjectBase<DerivedFA > & FA,
+      const Eigen::PlainObjectBase<DerivedVB > & VB,
+      const Eigen::PlainObjectBase<DerivedFB > & FB,
+      const MeshBooleanType & type,
+      Eigen::PlainObjectBase<DerivedVC > & VC,
+      Eigen::PlainObjectBase<DerivedFC > & FC,
+      Eigen::PlainObjectBase<DerivedJ > & J,
+      Eigen::PlainObjectBase<DerivedI > & I);
     template <
       typename DerivedVA,
-               typename DerivedFA,
-               typename DerivedVB,
-               typename DerivedFB,
-               typename DerivedVC,
-               typename DerivedFC,
-               typename DerivedJ>
-                 IGL_INLINE void mesh_boolean(
-                     const Eigen::PlainObjectBase<DerivedVA > & VA,
-                     const Eigen::PlainObjectBase<DerivedFA > & FA,
-                     const Eigen::PlainObjectBase<DerivedVB > & VB,
-                     const Eigen::PlainObjectBase<DerivedFB > & FB,
-                     const MeshBooleanType & type,
-                     const std::function<void(
-                       const Eigen::Matrix<
-                         typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-                       const Eigen::Matrix<
-                         typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
-                       Eigen::Matrix<
-                         typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-                       Eigen::Matrix<
-                         typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
-                       Eigen::Matrix<
-                         typename  DerivedJ::Scalar,Eigen::Dynamic,1>&)> 
-                     & resolve_fun,
-                     Eigen::PlainObjectBase<DerivedVC > & VC,
-                     Eigen::PlainObjectBase<DerivedFC > & FC,
-                     Eigen::PlainObjectBase<DerivedJ > & J);
+      typename DerivedFA,
+      typename DerivedVB,
+      typename DerivedFB,
+      typename DerivedVC,
+      typename DerivedFC>
+    IGL_INLINE void mesh_boolean(
+      const Eigen::PlainObjectBase<DerivedVA > & VA,
+      const Eigen::PlainObjectBase<DerivedFA > & FA,
+      const Eigen::PlainObjectBase<DerivedVB > & VB,
+      const Eigen::PlainObjectBase<DerivedFB > & FB,
+      const MeshBooleanType & type,
+      Eigen::PlainObjectBase<DerivedVC > & VC,
+      Eigen::PlainObjectBase<DerivedFC > & FC);
+
   }
 }
 

include/igl/boundary_facets.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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 "boundary_facets.h"
 #include "face_occurrences.h"
@@ -98,7 +98,6 @@ IGL_INLINE void igl::boundary_facets(
 
 }
 
-#ifndef IGL_NO_EIGEN
 #include "list_to_matrix.h"
 #include "matrix_to_list.h"
 
@@ -122,11 +121,10 @@ template <typename DerivedT, typename Ret>
 Ret igl::boundary_facets(
   const Eigen::PlainObjectBase<DerivedT>& T)
 {
-  Ret F; 
+  Ret F;
   igl::boundary_facets(T,F);
   return F;
 }
-#endif
 
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/boundary_facets.h

@@ -1,17 +1,15 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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_BOUNDARY_FACETS_H
 #define IGL_BOUNDARY_FACETS_H
 #include "igl_inline.h"
 
-#ifndef IGL_NO_EIGEN
-#  include <Eigen/Dense>
-#endif
+#include <Eigen/Dense>
 
 #include <vector>
 
@@ -34,7 +32,6 @@ namespace igl
     const std::vector<std::vector<IntegerT> > & T,
     std::vector<std::vector<IntegerF> > & F);
 
-#ifndef IGL_NO_EIGEN
   // Templates:
   //   DerivedT  integer-value: i.e. from MatrixXi
   //   DerivedF  integer-value: i.e. from MatrixXi
@@ -46,7 +43,6 @@ namespace igl
   template <typename DerivedT, typename Ret>
   Ret boundary_facets(
     const Eigen::PlainObjectBase<DerivedT>& T);
-#endif
 }
 
 #ifndef IGL_STATIC_LIBRARY
@@ -54,4 +50,3 @@ namespace igl
 #endif
 
 #endif
-

include/igl/cgal/is_inside.cpp

@@ -0,0 +1,374 @@
+#include "is_inside.h"
+
+#include <cassert>
+#include <list>
+#include <limits>
+#include <vector>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include "order_facets_around_edge.h"
+#include "assign_scalar.h"
+
+namespace igl {
+    namespace cgal {
+        namespace is_inside_helper {
+            typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+            typedef Kernel::Ray_3 Ray_3;
+            typedef Kernel::Point_3 Point_3;
+            typedef Kernel::Vector_3 Vector_3;
+            typedef Kernel::Triangle_3 Triangle;
+            typedef Kernel::Plane_3 Plane_3;
+            typedef std::vector<Triangle>::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;
+
+            enum ElementType { VERTEX, EDGE, FACE };
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            ElementType determine_element_type(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const size_t fid, const Point_3& p,
+                    size_t& element_index) {
+                const Eigen::Vector3i f = F.row(I(fid, 0));
+                const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+                const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+                const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+
+                if (p == p0) { element_index = 0; return VERTEX; }
+                if (p == p1) { element_index = 1; return VERTEX; }
+                if (p == p2) { element_index = 2; return VERTEX; }
+                if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+                if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+                if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+                element_index = 0;
+                return FACE;
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            void extract_adj_faces(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const size_t s, const size_t d,
+                    std::vector<int>& adj_faces) {
+                const size_t num_faces = I.rows();
+                for (size_t i=0; i<num_faces; i++) {
+                    Eigen::Vector3i f = F.row(I(i, 0));
+                    if ((f[0] == s && f[1] == d) ||
+                        (f[1] == s && f[2] == d) ||
+                        (f[2] == s && f[0] == d)) {
+                        adj_faces.push_back((I(i, 0)+1) * -1);
+                        continue;
+                    }
+                    if ((f[0] == d && f[1] == s) ||
+                        (f[1] == d && f[2] == s) ||
+                        (f[2] == d && f[0] == s)) {
+                        adj_faces.push_back(I(i, 0)+1);
+                        continue;
+                    }
+                }
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            void extract_adj_vertices(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const size_t v, std::vector<int>& adj_vertices) {
+                std::set<size_t> unique_adj_vertices;
+                const size_t num_faces = I.rows();
+                for (size_t i=0; i<num_faces; i++) {
+                    Eigen::Vector3i f = F.row(I(i, 0));
+                    assert((f.array() < V.rows()).all());
+                    if (f[0] == v) {
+                        unique_adj_vertices.insert(f[1]);
+                        unique_adj_vertices.insert(f[2]);
+                    } else if (f[1] == v) {
+                        unique_adj_vertices.insert(f[0]);
+                        unique_adj_vertices.insert(f[2]);
+                    } else if (f[2] == v) {
+                        unique_adj_vertices.insert(f[0]);
+                        unique_adj_vertices.insert(f[1]);
+                    }
+                }
+                adj_vertices.resize(unique_adj_vertices.size());
+                std::copy(unique_adj_vertices.begin(),
+                        unique_adj_vertices.end(),
+                        adj_vertices.begin());
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_edge_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t s, size_t d) {
+                // Algorithm:
+                //
+                // Order the adj faces around the edge (s,d) clockwise using
+                // query point as pivot.  (i.e. The first face of the ordering
+                // is directly after the pivot point, and the last face is
+                // directly before the pivot.)
+                //
+                // The point is outside if the first and last faces of the
+                // ordering forms a convex angle.  This check can be done
+                // without any construction by looking at the orientation of the
+                // faces.  The angle is convex iff the first face contains (s,d)
+                // as an edge and the last face contains (d,s) as an edge.
+                //
+                // The point is inside if the first and last faces of the
+                // ordering forms a concave angle.  That is the first face
+                // contains (d,s) as an edge and the last face contains (s,d) as
+                // an edge.
+                //
+                // In the special case of duplicated faces. I.e. multiple faces
+                // sharing the same 3 corners, but not necessarily the same
+                // orientation.  The ordering will always rank faces containing
+                // edge (s,d) before faces containing edge (d,s).
+                //
+                // Therefore, if there are any duplicates of the first faces,
+                // the ordering will always choose the one with edge (s,d) if
+                // possible.  The same for the last face.
+                //
+                // In the very degenerated case where the first and last face
+                // are duplicates, but with different orientations, it is
+                // equally valid to think the angle formed by them is either 0
+                // or 360 degrees.  By default, 0 degree is used, and thus the
+                // query point is outside.
+
+                std::vector<int> adj_faces;
+                extract_adj_faces(V, F, I, s, d, adj_faces);
+                const size_t num_adj_faces = adj_faces.size();
+                assert(num_adj_faces > 0);
+
+                DerivedV pivot_point(1, 3);
+                igl::cgal::assign_scalar(query.x(), pivot_point(0, 0));
+                igl::cgal::assign_scalar(query.y(), pivot_point(0, 1));
+                igl::cgal::assign_scalar(query.z(), pivot_point(0, 2));
+                Eigen::VectorXi order;
+                order_facets_around_edge(V, F, s, d,
+                        adj_faces, pivot_point, order);
+                assert(order.size() == num_adj_faces);
+                if (adj_faces[order[0]] > 0 &&
+                    adj_faces[order[num_adj_faces-1] < 0]) {
+                    return true;
+                } else if (adj_faces[order[0]] < 0 &&
+                    adj_faces[order[num_adj_faces-1] > 0]) {
+                    return false;
+                } else {
+                    throw "The input mesh does not represent a valid volume";
+                }
+                assert(false);
+                return false;
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_vertex_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t s) {
+                std::vector<int> adj_vertices;
+                extract_adj_vertices(V, F, I, s, adj_vertices);
+                const size_t num_adj_vertices = adj_vertices.size();
+
+                std::vector<Point_3> adj_points;
+                for (size_t i=0; i<num_adj_vertices; i++) {
+                    const size_t vi = adj_vertices[i];
+                    adj_points.emplace_back(V(vi,0), V(vi,1), V(vi,2));
+                }
+
+                // A plane is on the exterior if all adj_points lies on or to
+                // one side of the plane.
+                auto is_on_exterior = [&](const Plane_3& separator) {
+                    size_t positive=0;
+                    size_t negative=0;
+                    size_t coplanar=0;
+                    for (const auto& point : adj_points) {
+                        switch(separator.oriented_side(point)) {
+                            case CGAL::ON_POSITIVE_SIDE:
+                                positive++;
+                                break;
+                            case CGAL::ON_NEGATIVE_SIDE:
+                                negative++;
+                                break;
+                            case CGAL::ON_ORIENTED_BOUNDARY:
+                                coplanar++;
+                                break;
+                            default:
+                                assert(false);
+                        }
+                    }
+                    auto query_orientation = separator.oriented_side(query);
+                    bool r =
+                        (positive == 0 && query_orientation == CGAL::POSITIVE)
+                        ||
+                        (negative == 0 && query_orientation == CGAL::NEGATIVE);
+                    return r;
+                };
+
+                size_t d = std::numeric_limits<size_t>::max();
+                Point_3 p(V(s,0), V(s,1), V(s,2));
+                for (size_t i=0; i<num_adj_vertices; i++) {
+                    const size_t vi = adj_vertices[i];
+                    for (size_t j=i+1; j<num_adj_vertices; j++) {
+                        const size_t vj = adj_vertices[j];
+                        Plane_3 separator(p, adj_points[i], adj_points[j]);
+                        assert(!separator.is_degenerate());
+                        if (is_on_exterior(separator)) {
+                            d = vi;
+                            assert(!CGAL::collinear(p, adj_points[i], query));
+                            break;
+                        }
+                    }
+                    if (d < V.rows()) break;
+                }
+                if (d > V.rows()) {
+                    // All adj faces are coplanar, use the first edge.
+                    d = adj_vertices[0];
+                }
+                return determine_point_edge_orientation(V, F, I, query, s, d);
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_face_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t fid) {
+                // Algorithm: A point is on the inside of a face if the
+                // tetrahedron formed by them is negatively oriented.
+
+                Eigen::Vector3i f = F.row(I(fid, 0));
+                const Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+                const Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+                const Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+                auto result = CGAL::orientation(v0, v1, v2, query);
+                assert(result != CGAL::COPLANAR);
+                return result == CGAL::NEGATIVE;
+            }
+        }
+    }
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedI>
+IGL_INLINE bool igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V1,
+        const Eigen::PlainObjectBase<DerivedF>& F1,
+        const Eigen::PlainObjectBase<DerivedI>& I1,
+        const Eigen::PlainObjectBase<DerivedV>& V2,
+        const Eigen::PlainObjectBase<DerivedF>& F2,
+        const Eigen::PlainObjectBase<DerivedI>& I2) {
+    assert(F1.rows() > 0);
+    assert(I1.rows() > 0);
+    assert(F2.rows() > 0);
+    assert(I2.rows() > 0);
+
+    const Eigen::Vector3i& f = F1.row(I1(0, 0));
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> query(
+            (V1(f[0],0) + V1(f[1],0) + V1(f[2],0))/3.0,
+            (V1(f[0],1) + V1(f[1],1) + V1(f[2],1))/3.0,
+            (V1(f[0],2) + V1(f[1],2) + V1(f[2],2))/3.0);
+    Eigen::VectorXi inside;
+    igl::cgal::is_inside(V2, F2, I2, query, inside);
+    assert(inside.size() == 1);
+    return inside[0];
+}
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V1,
+        const Eigen::PlainObjectBase<DerivedF>& F1,
+        const Eigen::PlainObjectBase<DerivedV>& V2,
+        const Eigen::PlainObjectBase<DerivedF>& F2) {
+    Eigen::VectorXi I1(F1.rows()), I2(F2.rows());
+    I1.setLinSpaced(F1.rows(), 0, F1.rows()-1);
+    I2.setLinSpaced(F2.rows(), 0, F2.rows()-1);
+    return igl::cgal::is_inside(V1, F1, I1, V2, F2, I2);
+}
+
+template<typename DerivedV, typename DerivedF, typename DerivedI,
+    typename DerivedP, typename DerivedB>
+IGL_INLINE void igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside) {
+    using namespace igl::cgal::is_inside_helper;
+    assert(F.rows() > 0);
+    assert(I.rows() > 0);
+
+    const size_t num_faces = I.rows();
+    std::vector<Triangle> triangles;
+    for (size_t i=0; i<num_faces; i++) {
+        const Eigen::Vector3i f = F.row(I(i, 0));
+        triangles.emplace_back(
+                Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+                Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+                Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+        assert(!triangles.back().is_degenerate());
+    }
+    Tree tree(triangles.begin(), triangles.end());
+    tree.accelerate_distance_queries();
+
+    const size_t num_queries = P.rows();
+    inside.resize(num_queries);
+    for (size_t i=0; i<num_queries; i++) {
+        const Point_3 query(P(i,0), P(i,1), P(i,2));
+        auto projection = tree.closest_point_and_primitive(query);
+        auto closest_point = projection.first;
+        size_t fid = projection.second - triangles.begin();
+
+        size_t element_index;
+        switch (determine_element_type(
+                    V, F, I, fid, closest_point, element_index)) {
+            case VERTEX:
+                {
+                    const size_t s = F(I(fid, 0), element_index);
+                    inside[i] = determine_point_vertex_orientation(
+                            V, F, I, query, s);
+                }
+                break;
+            case EDGE:
+                {
+                    const size_t s = F(I(fid, 0), (element_index+1)%3);
+                    const size_t d = F(I(fid, 0), (element_index+2)%3);
+                    inside[i] = determine_point_edge_orientation(
+                            V, F, I, query, s, d);
+                }
+                break;
+            case FACE:
+                inside[i] = determine_point_face_orientation(V, F, I, query, fid);
+                break;
+            default:
+                assert(false);
+        }
+    }
+}
+
+template<typename DerivedV, typename DerivedF, typename DerivedP,
+    typename DerivedB>
+IGL_INLINE void igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside) {
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+    igl::cgal::is_inside(V, F, I, P, inside);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::cgal::is_inside<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::cgal::is_inside<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::cgal::is_inside<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/cgal/is_inside.h

@@ -0,0 +1,82 @@
+#ifndef IS_INSIDE
+#define IS_INSIDE
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+    namespace cgal {
+
+        // Determine if mesh (V1, F1) is inside of mesh (V2, F2).
+        //
+        // Precondition:
+        // Both mesh must represent closed, self-intersection free,
+        // non-degenerated surfaces that
+        // are the boundary of 3D volumes. In addition, (V1, F1) must not
+        // intersect with (V2, F2).
+        //
+        // Inputs:
+        //   V1  #V1 by 3 list of vertex position of mesh 1
+        //   F1  #F1 by 3 list of triangles indices into V1
+        //   I1  #I1 list of indices into F1, faces to consider.
+        //   V2  #V2 by 3 list of vertex position of mesh 2
+        //   F2  #F2 by 3 list of triangles indices into V2
+        //   I2  #I2 list of indices into F2, faces to consider.
+        //
+        // Outputs:
+        //   return true iff (V1, F1) is entirely inside of (V2, F2).
+        template<typename DerivedV, typename DerivedF, typename DerivedI>
+            IGL_INLINE bool is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V1,
+                    const Eigen::PlainObjectBase<DerivedF>& F1,
+                    const Eigen::PlainObjectBase<DerivedI>& I1,
+                    const Eigen::PlainObjectBase<DerivedV>& V2,
+                    const Eigen::PlainObjectBase<DerivedF>& F2,
+                    const Eigen::PlainObjectBase<DerivedI>& I2);
+
+        template<typename DerivedV, typename DerivedF>
+            IGL_INLINE bool is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V1,
+                    const Eigen::PlainObjectBase<DerivedF>& F1,
+                    const Eigen::PlainObjectBase<DerivedV>& V2,
+                    const Eigen::PlainObjectBase<DerivedF>& F2);
+
+        // Determine if queries points are inside of mesh (V, F).
+        //
+        // Precondition:
+        // The input mesh must be a closed, self-intersection free,
+        // non-degenerated surface.  Queries points must be either inside or
+        // outside of the mesh.
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertex positions.
+        //   F  #F by 3 array of triangles.
+        //   I  #I list of triangle indices to consider.
+        //   P  #P by 3 array of query points.
+        //
+        // Outputs:
+        //   inside  #P list of booleans that is true iff the corresponding
+        //           query point is inside of the mesh.
+        template<typename DerivedV, typename DerivedF, typename DerivedI,
+            typename DerivedP, typename DerivedB>
+            IGL_INLINE void is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Eigen::PlainObjectBase<DerivedP>& P,
+                    Eigen::PlainObjectBase<DerivedB>& inside);
+
+        template<typename DerivedV, typename DerivedF, typename DerivedP,
+            typename DerivedB>
+            IGL_INLINE void is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedP>& P,
+                    Eigen::PlainObjectBase<DerivedB>& inside);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "is_inside.cpp"
+#endif
+#endif

include/igl/cgal/order_facets_around_edge.cpp

@@ -1,16 +1,23 @@
 #include "order_facets_around_edge.h"
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
-namespace igl {
-    namespace cgal {
-        namespace order_facets_around_edges_helper {
+namespace igl
+{
+    namespace cgal
+    {
+        namespace order_facets_around_edges_helper
+        {
             template<typename T>
-            std::vector<size_t> index_sort(const std::vector<T>& data) {
+            std::vector<size_t> index_sort(const std::vector<T>& data)
+            {
                 const size_t len = data.size();
                 std::vector<size_t> order(len);
-                for (size_t i=0; i<len; i++) order[i] = i;
-
-                auto comp = [&](size_t i, size_t j) {
+                for (size_t i=0; i<len; i++) 
+                {
+                    order[i] = i;
+                }
+                auto comp = [&](size_t i, size_t j) 
+                {
                     return data[i] < data[j];
                 };
                 std::sort(order.begin(), order.end(), comp);
@@ -28,10 +35,11 @@ template<
 void igl::cgal::order_facets_around_edge(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
-    size_t s, size_t d, 
+    size_t s,
+    size_t d, 
     const std::vector<int>& adj_faces,
-    Eigen::PlainObjectBase<DerivedI>& order) {
-
+    Eigen::PlainObjectBase<DerivedI>& order)
+{
     using namespace igl::cgal::order_facets_around_edges_helper;
 
     // Although we only need exact predicates in the algorithm,
@@ -41,37 +49,42 @@ void igl::cgal::order_facets_around_edge(
     typedef K::Point_3 Point_3;
     typedef K::Plane_3 Plane_3;
 
-    auto get_face_index = [&](int adj_f)->size_t{
+    auto get_face_index = [&](int adj_f)->size_t
+    {
         return abs(adj_f) - 1;
     };
 
-    auto get_opposite_vertex = [&](size_t fid)->size_t {
-        if (F(fid, 0) != s && F(fid, 0) != d) return F(fid, 0);
-        if (F(fid, 1) != s && F(fid, 1) != d) return F(fid, 1);
-        if (F(fid, 2) != s && F(fid, 2) != d) return F(fid, 2);
+    auto get_opposite_vertex = [&](size_t fid)->size_t
+    {
+        typedef typename DerivedF::Scalar Index;
+        if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
+        if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
+        if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
         assert(false);
         return -1;
     };
 
     // Handle base cases
-    if (adj_faces.size() == 0) {
+    if (adj_faces.size() == 0) 
+    {
         order.resize(0, 1);
         return;
-    } else if (adj_faces.size() == 1) {
+    } else if (adj_faces.size() == 1)
+    {
         order.resize(1, 1);
         order(0, 0) = 0;
         return;
-    } else if (adj_faces.size() == 2) {
-        const size_t o1 =
-            get_opposite_vertex(get_face_index(adj_faces[0]));
-        const size_t o2 =
-            get_opposite_vertex(get_face_index(adj_faces[1]));
+    } else if (adj_faces.size() == 2)
+    {
+        const size_t o1 = get_opposite_vertex(get_face_index(adj_faces[0]));
+        const size_t o2 = get_opposite_vertex(get_face_index(adj_faces[1]));
         const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
         const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
         const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
         const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
         order.resize(2, 1);
-        switch (CGAL::orientation(ps, pd, p1, p2)) {
+        switch (CGAL::orientation(ps, pd, p1, p2))
+        {
             case CGAL::POSITIVE:
                 order(0, 0) = 1;
                 order(1, 0) = 0;
@@ -82,19 +95,27 @@ void igl::cgal::order_facets_around_edge(
                 break;
             case CGAL::COPLANAR:
                 {
-                    Plane_3 P1(ps, pd, p1);
-                    Plane_3 P2(ps, pd, p2);
-                    if (P1.orthogonal_direction() == P2.orthogonal_direction()){
-                        // Duplicated face, use index to break tie.
-                        order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
-                        order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
-                    } else {
-                        // Coplanar faces, one on each side of the edge.
-                        // It is equally valid to order them (0, 1) or (1, 0).
-                        // I cannot think of any reason to prefer one to the
-                        // other.  So just use (0, 1) ordering by default.
-                        order(0, 0) = 0;
-                        order(1, 0) = 1;
+                    switch (CGAL::coplanar_orientation(ps, pd, p1, p2)) {
+                        case CGAL::POSITIVE:
+                            // Duplicated face, use index to break tie.
+                            order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
+                            order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+                            break;
+                        case CGAL::NEGATIVE:
+                            // Coplanar faces, one on each side of the edge.
+                            // It is equally valid to order them (0, 1) or (1, 0).
+                            // I cannot think of any reason to prefer one to the
+                            // other.  So just use (0, 1) ordering by default.
+                            order(0, 0) = 0;
+                            order(1, 0) = 1;
+                            break;
+                        case CGAL::COLLINEAR:
+                            std::cerr << "Degenerated triangle detected." <<
+                                std::endl;
+                            assert(false);
+                            break;
+                        default:
+                            assert(false);
                     }
                 }
                 break;
@@ -105,8 +126,7 @@ void igl::cgal::order_facets_around_edge(
     }
 
     const size_t num_adj_faces = adj_faces.size();
-    const size_t o = get_opposite_vertex(
-            get_face_index(adj_faces[0]));
+    const size_t o = get_opposite_vertex( get_face_index(adj_faces[0]));
     const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
     const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
     const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
@@ -114,9 +134,9 @@ void igl::cgal::order_facets_around_edge(
     assert(!separator.is_degenerate());
 
     std::vector<Point_3> opposite_vertices;
-    for (size_t i=0; i<num_adj_faces; i++) {
-        const size_t o = get_opposite_vertex(
-                get_face_index(adj_faces[i]));
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
+        const size_t o = get_opposite_vertex( get_face_index(adj_faces[i]));
         opposite_vertices.emplace_back(
                 V(o, 0), V(o, 1), V(o, 2));
     }
@@ -131,7 +151,8 @@ void igl::cgal::order_facets_around_edge(
     std::vector<size_t> tie_positive_oriented_index;
     std::vector<size_t> tie_negative_oriented_index;
 
-    for (size_t i=0; i<num_adj_faces; i++) {
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
         const int f = adj_faces[i];
         const Point_3& p_a = opposite_vertices[i];
         auto orientation = separator.oriented_side(p_a);
@@ -146,17 +167,21 @@ void igl::cgal::order_facets_around_edge(
                 break;
             case CGAL::ON_ORIENTED_BOUNDARY:
                 {
-                    const Plane_3 other(p_s, p_d, p_a);
-                    const auto target_dir = separator.orthogonal_direction();
-                    const auto query_dir = other.orthogonal_direction();
-                    if (target_dir == query_dir) {
-                        tie_positive_oriented.push_back(f);
-                        tie_positive_oriented_index.push_back(i);
-                    } else if (target_dir == -query_dir) {
-                        tie_negative_oriented.push_back(f);
-                        tie_negative_oriented_index.push_back(i);
-                    } else {
-                        assert(false);
+                    auto inplane_orientation = CGAL::coplanar_orientation(
+                            p_s, p_d, p_o, p_a);
+                    switch (inplane_orientation) {
+                        case CGAL::POSITIVE:
+                            tie_positive_oriented.push_back(f);
+                            tie_positive_oriented_index.push_back(i);
+                            break;
+                        case CGAL::NEGATIVE:
+                            tie_negative_oriented.push_back(f);
+                            tie_negative_oriented_index.push_back(i);
+                            break;
+                        case CGAL::COLLINEAR:
+                        default:
+                            assert(false);
+                            break;
                     }
                 }
                 break;
@@ -169,10 +194,8 @@ void igl::cgal::order_facets_around_edge(
     Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
     order_facets_around_edge(V, F, s, d, positive_side, positive_order);
     order_facets_around_edge(V, F, s, d, negative_side, negative_order);
-    std::vector<size_t> tie_positive_order =
-        index_sort(tie_positive_oriented);
-    std::vector<size_t> tie_negative_order =
-        index_sort(tie_negative_oriented);
+    std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
+    std::vector<size_t> tie_negative_order = index_sort(tie_negative_oriented);
 
     // Copy results into order vector.
     const size_t tie_positive_size = tie_positive_oriented.size();
@@ -180,28 +203,30 @@ void igl::cgal::order_facets_around_edge(
     const size_t positive_size = positive_order.size();
     const size_t negative_size = negative_order.size();
 
-    order.resize(tie_positive_size + positive_size +
-            tie_negative_size + negative_size, 1);
+    order.resize(
+      tie_positive_size + positive_size + tie_negative_size + negative_size,1);
 
     size_t count=0;
-    for (size_t i=0; i<tie_positive_size; i++) {
-        order(count+i, 0) =
-            tie_positive_oriented_index[tie_positive_order[i]];
+    for (size_t i=0; i<tie_positive_size; i++)
+    {
+        order(count+i, 0) = tie_positive_oriented_index[tie_positive_order[i]];
     }
     count += tie_positive_size;
 
-    for (size_t i=0; i<negative_size; i++) {
+    for (size_t i=0; i<negative_size; i++) 
+    {
         order(count+i, 0) = negative_side_index[negative_order(i, 0)];
     }
     count += negative_size;
 
-    for (size_t i=0; i<tie_negative_size; i++) {
-        order(count+i, 0) =
-            tie_negative_oriented_index[tie_negative_order[i]];
+    for (size_t i=0; i<tie_negative_size; i++)
+    {
+        order(count+i, 0) = tie_negative_oriented_index[tie_negative_order[i]];
     }
     count += tie_negative_size;
 
-    for (size_t i=0; i<positive_size; i++) {
+    for (size_t i=0; i<positive_size; i++)
+    {
         order(count+i, 0) = positive_side_index[positive_order(i, 0)];
     }
     count += positive_size;
@@ -209,19 +234,22 @@ void igl::cgal::order_facets_around_edge(
 
     // Find the correct start point.
     size_t start_idx = 0;
-    for (size_t i=0; i<num_adj_faces; i++) {
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
         const Point_3& p_a = opposite_vertices[order(i, 0)];
         const Point_3& p_b =
             opposite_vertices[order((i+1)%num_adj_faces, 0)];
         auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
-        if (orientation == CGAL::POSITIVE) {
+        if (orientation == CGAL::POSITIVE)
+        {
             // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
             // more than 180 degrees.
             start_idx = (i+1)%num_adj_faces;
             break;
         } else if (orientation == CGAL::COPLANAR &&
                 Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
-                Plane_3(p_s, p_d, p_b).orthogonal_direction()) {
+                Plane_3(p_s, p_d, p_b).orthogonal_direction())
+        {
             // All 4 points are coplanar, but p_a and p_b are on each side of
             // the edge (p_s, p_d).  This means the angle between triangle
             // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
@@ -230,7 +258,8 @@ void igl::cgal::order_facets_around_edge(
         }
     }
     DerivedI circular_order = order;
-    for (size_t i=0; i<num_adj_faces; i++) {
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
         order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
     }
 }
@@ -243,31 +272,55 @@ IGL_INLINE
 void igl::cgal::order_facets_around_edge(
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
-        size_t s, size_t d, 
+        size_t s,
+        size_t d, 
         const std::vector<int>& adj_faces,
         const Eigen::PlainObjectBase<DerivedV>& pivot_point,
-        Eigen::PlainObjectBase<DerivedI>& order) {
+        Eigen::PlainObjectBase<DerivedI>& order)
+{
 
     assert(V.cols() == 3);
     assert(F.cols() == 3);
-    auto signed_index_to_index = [&](int signed_idx) {
+    auto signed_index_to_index = [&](int signed_idx)
+    {
         return abs(signed_idx) -1;
     };
-    auto get_opposite_vertex_index = [&](size_t fid) {
-        if (F(fid, 0) != s && F(fid, 0) != d) return F(fid, 0);
-        if (F(fid, 1) != s && F(fid, 1) != d) return F(fid, 1);
-        if (F(fid, 2) != s && F(fid, 2) != d) return F(fid, 2);
+    auto get_opposite_vertex_index = [&](size_t fid)
+    {
+        typedef typename DerivedF::Scalar Index;
+        if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
+        if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
+        if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
         assert(false);
         // avoid warning
         return -1;
     };
 
+    {
+        // Check if s, d and pivot are collinear.
+        typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+        K::Point_3 ps(V(s,0), V(s,1), V(s,2));
+        K::Point_3 pd(V(d,0), V(d,1), V(d,2));
+        K::Point_3 pp(pivot_point(0,0), pivot_point(0,1), pivot_point(0,2));
+        if (CGAL::collinear(ps, pd, pp)) {
+            throw "Pivot point is collinear with the outer edge!";
+        }
+    }
+
     const size_t N = adj_faces.size();
     const size_t num_faces = N + 1; // N adj faces + 1 pivot face
 
+    // Because face indices are used for tie breaking, the original face indices
+    // in the new faces array must be ascending.
+    auto comp = [&](int i, int j) {
+        return signed_index_to_index(i) < signed_index_to_index(j);
+    };
+    std::vector<int> ordered_adj_faces(adj_faces);
+    std::sort(ordered_adj_faces.begin(), ordered_adj_faces.end(), comp);
+
     DerivedV vertices(num_faces + 2, 3);
     for (size_t i=0; i<N; i++) {
-        const size_t fid = signed_index_to_index(adj_faces[i]);
+        const size_t fid = signed_index_to_index(ordered_adj_faces[i]);
         vertices.row(i) = V.row(get_opposite_vertex_index(fid));
     }
     vertices.row(N  ) = pivot_point;
@@ -275,8 +328,9 @@ void igl::cgal::order_facets_around_edge(
     vertices.row(N+2) = V.row(d);
 
     DerivedF faces(num_faces, 3);
-    for (size_t i=0; i<N; i++) {
-        if (adj_faces[i] < 0) {
+    for (size_t i=0; i<N; i++)
+    {
+        if (ordered_adj_faces[i] < 0) {
             faces(i,0) = N+1; // s
             faces(i,1) = N+2; // d
             faces(i,2) = i  ;
@@ -292,10 +346,13 @@ void igl::cgal::order_facets_around_edge(
     faces(N, 2) = N;
 
     std::vector<int> adj_faces_with_pivot(num_faces);
-    for (size_t i=0; i<num_faces; i++) {
-        if (faces(i,0) == N+1 && faces(i,1) == N+2) {
+    for (size_t i=0; i<num_faces; i++)
+    {
+        if (faces(i,0) == N+1 && faces(i,1) == N+2)
+        {
             adj_faces_with_pivot[i] = int(i+1) * -1;
-        } else {
+        } else
+        {
             adj_faces_with_pivot[i] = int(i+1);
         }
     }
@@ -308,15 +365,18 @@ void igl::cgal::order_facets_around_edge(
     assert(order_with_pivot.size() == num_faces);
     order.resize(N);
     size_t pivot_index = num_faces + 1;
-    for (size_t i=0; i<num_faces; i++) {
-        if (order_with_pivot[i] == N) {
+    for (size_t i=0; i<num_faces; i++)
+    {
+        if (order_with_pivot[i] == N)
+        {
             pivot_index = i;
             break;
         }
     }
     assert(pivot_index < num_faces);
 
-    for (size_t i=0; i<N; i++) {
+    for (size_t i=0; i<N; i++)
+    {
         order[i] = order_with_pivot[(pivot_index+i+1)%num_faces];
     }
 }

include/igl/cgal/order_facets_around_edges.cpp

@@ -15,9 +15,7 @@ template<
     typename DerivedV,
     typename DerivedF,
     typename DerivedN,
-    typename DerivedE,
     typename DeriveduE,
-    typename DerivedEMAP,
     typename uE2EType,
     typename uE2oEType,
     typename uE2CType >
@@ -28,9 +26,7 @@ igl::cgal::order_facets_around_edges(
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedN>& N,
-        const Eigen::PlainObjectBase<DerivedE>& E,
         const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
         const std::vector<std::vector<uE2EType> >& uE2E,
         std::vector<std::vector<uE2oEType> >& uE2oE,
         std::vector<std::vector<uE2CType > >& uE2C ) {
@@ -159,9 +155,7 @@ template<
     typename DerivedV,
     typename DerivedF,
     typename DerivedN,
-    typename DerivedE,
     typename DeriveduE,
-    typename DerivedEMAP,
     typename uE2EType,
     typename uE2oEType,
     typename uE2CType >
@@ -172,9 +166,7 @@ igl::cgal::order_facets_around_edges(
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedN>& N,
-        const Eigen::PlainObjectBase<DerivedE>& E,
         const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
         const std::vector<std::vector<uE2EType> >& uE2E,
         std::vector<std::vector<uE2oEType> >& uE2oE,
         std::vector<std::vector<uE2CType > >& uE2C ) {
@@ -257,23 +249,19 @@ igl::cgal::order_facets_around_edges(
 template<
     typename DerivedV,
     typename DerivedF,
-    typename DerivedE,
     typename DeriveduE,
-    typename DerivedEMAP,
     typename uE2EType,
     typename uE2oEType,
     typename uE2CType >
 IGL_INLINE void igl::cgal::order_facets_around_edges(
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedE>& E,
         const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
         const std::vector<std::vector<uE2EType> >& uE2E,
         std::vector<std::vector<uE2oEType> >& uE2oE,
         std::vector<std::vector<uE2CType > >& uE2C ) {
 
-    typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
+    //typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
     const size_t num_faces = F.rows();
     const size_t num_undirected_edges = uE.rows();
 
@@ -296,7 +284,7 @@ IGL_INLINE void igl::cgal::order_facets_around_edges(
         const auto ref_corner_s = (ref_corner_o+1)%3;
         const auto ref_corner_d = (ref_corner_o+2)%3;
 
-        const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
+        //const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
         const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
         const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
 
@@ -332,7 +320,7 @@ IGL_INLINE void igl::cgal::order_facets_around_edges(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
-template void igl::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, bool>(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, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
-template void igl::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, bool>(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, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
+template void igl::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
+template void igl::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(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, 2, 0, -1, 2> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
+template void igl::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(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, 2, 0, -1, 2> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
 #endif

include/igl/cgal/order_facets_around_edges.h

@@ -29,9 +29,7 @@ namespace igl
     //   V    #V by 3 list of vertices.
     //   F    #F by 3 list of faces
     //   N    #F by 3 list of face normals.
-    //   E    #F*3 by 2 list vertex indices, represents directed edges.
     //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
-    //  EMAP  #F*3 list of indices that maps E to uE. (a many-to-one map)
     //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
     //
     // Outputs:
@@ -44,9 +42,7 @@ namespace igl
         typename DerivedV,
         typename DerivedF,
         typename DerivedN,
-        typename DerivedE,
         typename DeriveduE,
-        typename DerivedEMAP,
         typename uE2EType,
         typename uE2oEType,
         typename uE2CType >
@@ -57,9 +53,7 @@ namespace igl
             const Eigen::PlainObjectBase<DerivedV>& V,
             const Eigen::PlainObjectBase<DerivedF>& F,
             const Eigen::PlainObjectBase<DerivedN>& N,
-            const Eigen::PlainObjectBase<DerivedE>& E,
             const Eigen::PlainObjectBase<DeriveduE>& uE,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
             const std::vector<std::vector<uE2EType> >& uE2E,
             std::vector<std::vector<uE2oEType> >& uE2oE,
             std::vector<std::vector<uE2CType > >& uE2C );
@@ -68,9 +62,7 @@ namespace igl
         typename DerivedV,
         typename DerivedF,
         typename DerivedN,
-        typename DerivedE,
         typename DeriveduE,
-        typename DerivedEMAP,
         typename uE2EType,
         typename uE2oEType,
         typename uE2CType >
@@ -81,9 +73,7 @@ namespace igl
             const Eigen::PlainObjectBase<DerivedV>& V,
             const Eigen::PlainObjectBase<DerivedF>& F,
             const Eigen::PlainObjectBase<DerivedN>& N,
-            const Eigen::PlainObjectBase<DerivedE>& E,
             const Eigen::PlainObjectBase<DeriveduE>& uE,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
             const std::vector<std::vector<uE2EType> >& uE2E,
             std::vector<std::vector<uE2oEType> >& uE2oE,
             std::vector<std::vector<uE2CType > >& uE2C );
@@ -93,18 +83,14 @@ namespace igl
     template<
         typename DerivedV,
         typename DerivedF,
-        typename DerivedE,
         typename DeriveduE,
-        typename DerivedEMAP,
         typename uE2EType,
         typename uE2oEType,
         typename uE2CType >
     IGL_INLINE void order_facets_around_edges(
             const Eigen::PlainObjectBase<DerivedV>& V,
             const Eigen::PlainObjectBase<DerivedF>& F,
-            const Eigen::PlainObjectBase<DerivedE>& E,
             const Eigen::PlainObjectBase<DeriveduE>& uE,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
             const std::vector<std::vector<uE2EType> >& uE2E,
             std::vector<std::vector<uE2oEType> >& uE2oE,
             std::vector<std::vector<uE2CType > >& uE2C );

include/igl/cgal/outer_facet.cpp

@@ -41,9 +41,8 @@ IGL_INLINE void igl::cgal::outer_facet(
     //    If it has zero normal component, it is facing outside if it contains
     //    directed edge (s, d).  
 
-    typedef typename DerivedV::Scalar Scalar;
+    //typedef typename DerivedV::Scalar Scalar;
     typedef typename DerivedV::Index Index;
-    const size_t INVALID = std::numeric_limits<size_t>::max();
 
     Index s,d;
     Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;

include/igl/cgal/outer_hull.cpp

@@ -5,6 +5,7 @@
 // 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 "is_inside.h"
 #include "outer_hull.h"
 #include "order_facets_around_edges.h"
 #include "outer_facet.h"
@@ -17,6 +18,7 @@
 #include "../per_face_normals.h"
 #include "../writePLY.h"
 #include "../sort_angles.h"
+#include "../writePLY.h"
 
 #include <Eigen/Geometry>
 #include <vector>
@@ -48,7 +50,7 @@ IGL_INLINE void igl::cgal::outer_hull(
   typedef typename DerivedF::Index Index;
   Matrix<Index,DerivedF::RowsAtCompileTime,1> C;
   typedef Matrix<typename DerivedV::Scalar,Dynamic,DerivedV::ColsAtCompileTime> MatrixXV;
-  typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
+  //typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
   typedef Matrix<typename DerivedG::Scalar,Dynamic,DerivedG::ColsAtCompileTime> MatrixXG;
   typedef Matrix<typename DerivedJ::Scalar,Dynamic,DerivedJ::ColsAtCompileTime> MatrixXJ;
   const Index m = F.rows();
@@ -74,7 +76,7 @@ IGL_INLINE void igl::cgal::outer_hull(
 #endif
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
   typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
-  typedef Matrix<typename DerivedV::Scalar, 3, 1> Vector3F;
+  //typedef Matrix<typename DerivedV::Scalar, 3, 1> Vector3F;
   MatrixX2I E,uE;
   VectorXI EMAP;
   vector<vector<typename DerivedF::Index> > uE2E;
@@ -91,7 +93,7 @@ IGL_INLINE void igl::cgal::outer_hull(
 
   std::vector<std::vector<typename DerivedF::Index> > uE2oE;
   std::vector<std::vector<bool> > uE2C;
-  order_facets_around_edges(V, F, E, uE, EMAP, uE2E, uE2oE, uE2C);
+  order_facets_around_edges(V, F, uE, uE2E, uE2oE, uE2C);
   uE2E = uE2oE;
   VectorXI diIM(3*m);
   for (auto ue : uE2E) {
@@ -325,19 +327,19 @@ IGL_INLINE void igl::cgal::outer_hull(
 
   // Is A inside B? Assuming A and B are consistently oriented but closed and
   // non-intersecting.
-  const auto & is_component_inside_other = [](
-    const Eigen::MatrixXd & V,
+  const auto & has_overlapping_bbox = [](
+    const Eigen::PlainObjectBase<DerivedV> & V,
     const MatrixXV & BC,
     const MatrixXG & A,
     const MatrixXJ & AJ,
     const MatrixXG & B)->bool
   {
     const auto & bounding_box = [](
-      const Eigen::MatrixXd & V,
+      const Eigen::PlainObjectBase<DerivedV> & V,
       const MatrixXG & F)->
-      Eigen::MatrixXd
+        DerivedV
     {
-      Eigen::MatrixXd BB(2,3);
+      DerivedV BB(2,3);
       BB<<
          1e26,1e26,1e26,
         -1e26,-1e26,-1e26;
@@ -346,75 +348,73 @@ IGL_INLINE void igl::cgal::outer_hull(
       {
         for(size_t c = 0;c<3;c++)
         {
-          const auto & vfc = V.row(F(f,c));
-          BB.row(0) = BB.row(0).array().min(vfc.array()).eval();
-          BB.row(1) = BB.row(1).array().max(vfc.array()).eval();
+          const auto & vfc = V.row(F(f,c)).eval();
+          BB(0,0) = std::min(BB(0,0), vfc(0,0));
+          BB(0,1) = std::min(BB(0,1), vfc(0,1));
+          BB(0,2) = std::min(BB(0,2), vfc(0,2));
+          BB(1,0) = std::max(BB(1,0), vfc(0,0));
+          BB(1,1) = std::max(BB(1,1), vfc(0,1));
+          BB(1,2) = std::max(BB(1,2), vfc(0,2));
         }
       }
       return BB;
     };
     // A lot of the time we're dealing with unrelated, distant components: cull
     // them.
-    Eigen::MatrixXd ABB = bounding_box(V,A);
-    Eigen::MatrixXd BBB = bounding_box(V,B);
+    DerivedV ABB = bounding_box(V,A);
+    DerivedV BBB = bounding_box(V,B);
     if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
         (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
     {
       // bounding boxes do not overlap
       return false;
+    } else {
+      return true;
     }
-    ////////////////////////////////////////////////////////////////////////
-    // POTENTIAL ROBUSTNESS WEAK AREA
-    ////////////////////////////////////////////////////////////////////////
-    //
-
-    // winding_number_3 expects colmajor
-    // q could be so close (<~1e-15) to B that the winding number is not a robust way to
-    // determine inside/outsideness. We could try to find a _better_ q which is
-    // farther away, but couldn't they all be bad?
-    double q[3] = {
-        CGAL::to_double(BC(AJ(0), 0)),
-        CGAL::to_double(BC(AJ(0), 1)),
-        CGAL::to_double(BC(AJ(0), 2)) };
-    // In a perfect world, it's enough to test a single point.
-    double w;
-    winding_number_3(
-      V.data(),V.rows(),
-      B.data(),B.rows(),
-      q,1,&w);
-    return w > 0.5 || w < -0.5;
   };
 
-  Eigen::MatrixXd Vcol(V.rows(), V.cols());
-  for (size_t i=0; i<(size_t)V.rows(); i++) {
-      for (size_t j=0; j<(size_t)V.cols(); j++) {
-          Vcol(i, j) = CGAL::to_double(V(i, j));
-      }
-  }
-
   // Reject components which are completely inside other components
   vector<bool> keep(ncc,true);
   size_t nG = 0;
   // This is O( ncc * ncc * m)
   for(size_t id = 0;id<ncc;id++)
   {
+    if (!keep[id]) continue;
+    std::vector<size_t> unresolved;
     for(size_t oid = 0;oid<ncc;oid++)
     {
-      if(id == oid)
+      if(id == oid || !keep[oid])
       {
         continue;
       }
-      const bool inside = is_component_inside_other(Vcol,BC,vG[id],vJ[id],vG[oid]);
-#ifdef IGL_OUTER_HULL_DEBUG
-      cout<<id<<" is inside "<<oid<<" ? "<<inside<<endl;
-#endif
-      keep[id] = keep[id] && !inside;
+      if (has_overlapping_bbox(V, BC, vG[id], vJ[id], vG[oid])) {
+          unresolved.push_back(oid);
+      }
     }
-    if(keep[id])
-    {
-      nG += vJ[id].rows();
+    const size_t num_unresolved_components = unresolved.size();
+    DerivedV query_points(num_unresolved_components, 3);
+    for (size_t i=0; i<num_unresolved_components; i++) {
+        const size_t oid = unresolved[i];
+        DerivedF f = vG[oid].row(0);
+        query_points(i,0) = (V(f(0,0), 0) + V(f(0,1), 0) + V(f(0,2), 0))/3.0;
+        query_points(i,1) = (V(f(0,0), 1) + V(f(0,1), 1) + V(f(0,2), 1))/3.0;
+        query_points(i,2) = (V(f(0,0), 2) + V(f(0,1), 2) + V(f(0,2), 2))/3.0;
+    }
+    Eigen::VectorXi inside;
+    igl::cgal::is_inside(V, vG[id], query_points, inside);
+    assert(inside.size() == num_unresolved_components);
+    for (size_t i=0; i<num_unresolved_components; i++) {
+        if (inside[i]) {
+            const size_t oid = unresolved[i];
+            keep[oid] = false;
+        }
     }
   }
+  for (size_t id = 0; id<ncc; id++) {
+      if (keep[id]) {
+          nG += vJ[id].rows();
+      }
+  }
 
   // collect G and J across components
   G.resize(nG,3);

include/igl/cgal/projected_delaunay.cpp

@@ -32,10 +32,10 @@ IGL_INLINE void igl::cgal::projected_delaunay(
   typedef CGAL::Segment_3<Kernel>  Segment_3; 
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
   typedef CGAL::Plane_3<Kernel>    Plane_3;
-  typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
+  //typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
   typedef CGAL::Point_2<Kernel>    Point_2;
-  typedef CGAL::Segment_2<Kernel>  Segment_2; 
-  typedef CGAL::Triangle_2<Kernel> Triangle_2; 
+  //typedef CGAL::Segment_2<Kernel>  Segment_2; 
+  //typedef CGAL::Triangle_2<Kernel> Triangle_2; 
   typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
   typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
   typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;

include/igl/cgal/remesh_self_intersections.cpp

@@ -88,4 +88,5 @@ template void igl::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3,
 template void igl::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::cgal::remesh_self_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/conjugate_frame_fields.h

@@ -6,8 +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/.
 
-#ifndef IGL_CONJUGATE_FRAME_FIELDS
-#define IGL_CONJUGATE_FRAME_FIELDS
+#ifndef IGL_CONJUGATE_FRAME_FIELDS_H
+#define IGL_CONJUGATE_FRAME_FIELDS_H
 #include "igl_inline.h"
 #include "ConjugateFFSolverData.h"
 #include <Eigen/Core>
@@ -60,4 +60,4 @@ namespace igl {
 #endif
 
 
-#endif /* defined(IGL_CONJUGATE_FRAME_FIELDS) */
+#endif

include/igl/outer_element.cpp

@@ -23,58 +23,59 @@ IGL_INLINE void igl::outer_vertex(
         IndexType & v_index,
         Eigen::PlainObjectBase<DerivedA> & A)
 {
-  // Algorithm: 
-  //    Find an outer vertex (i.e. vertex reachable from infinity)
-  //    Return the vertex with the largest X value.
-  //    If there is a tie, pick the one with largest Y value.
-  //    If there is still a tie, pick the one with the largest Z value.
-  //    If there is still a tie, then there are duplicated vertices within the
-  //    mesh, which violates the precondition.
-  const size_t INVALID = std::numeric_limits<size_t>::max();
-  const size_t num_selected_faces = I.rows();
-  std::vector<size_t> candidate_faces;
-  size_t outer_vid = INVALID;
-  typename DerivedV::Scalar outer_val = 0;
-  for (size_t i=0; i<num_selected_faces; i++)
-  {
-    size_t f = I(i);
-    for (size_t j=0; j<3; j++)
+    // Algorithm: 
+    //    Find an outer vertex (i.e. vertex reachable from infinity)
+    //    Return the vertex with the largest X value.
+    //    If there is a tie, pick the one with largest Y value.
+    //    If there is still a tie, pick the one with the largest Z value.
+    //    If there is still a tie, then there are duplicated vertices within the
+    //    mesh, which violates the precondition.
+    typedef typename DerivedF::Scalar Index;
+    const Index INVALID = std::numeric_limits<Index>::max();
+    const size_t num_selected_faces = I.rows();
+    std::vector<size_t> candidate_faces;
+    Index outer_vid = INVALID;
+    typename DerivedV::Scalar outer_val = 0;
+    for (size_t i=0; i<num_selected_faces; i++)
     {
-      auto v = F(f, j);
-      auto vx = V(v, 0);
-      if (outer_vid == INVALID || vx > outer_val)
-      {
-        outer_val = vx;
-        outer_vid = v;
-        candidate_faces = {f};
-      } else if (v == outer_vid)
-      {
-        candidate_faces.push_back(f);
-      } else if (vx == outer_val)
-      {
-        // Break tie.
-        auto vy = V(v,1);
-        auto vz = V(v, 2);
-        auto outer_y = V(outer_vid, 1);
-        auto outer_z = V(outer_vid, 2);
-        assert(!(vy == outer_y && vz == outer_z));
-        bool replace = (vy > outer_y) ||
-          ((vy == outer_y) && (vz > outer_z));
-        if (replace)
+        size_t f = I(i);
+        for (size_t j=0; j<3; j++)
         {
-          outer_val = vx;
-          outer_vid = v;
-          candidate_faces = {f};
+            Index v = F(f, j);
+            auto vx = V(v, 0);
+            if (outer_vid == INVALID || vx > outer_val)
+            {
+                outer_val = vx;
+                outer_vid = v;
+                candidate_faces = {f};
+            } else if (v == outer_vid)
+            {
+                candidate_faces.push_back(f);
+            } else if (vx == outer_val)
+            {
+                // Break tie.
+                auto vy = V(v,1);
+                auto vz = V(v, 2);
+                auto outer_y = V(outer_vid, 1);
+                auto outer_z = V(outer_vid, 2);
+                assert(!(vy == outer_y && vz == outer_z));
+                bool replace = (vy > outer_y) ||
+                    ((vy == outer_y) && (vz > outer_z));
+                if (replace)
+                {
+                    outer_val = vx;
+                    outer_vid = v;
+                    candidate_faces = {f};
+                }
+            }
         }
-      }
     }
-  }
 
-  assert(outer_vid != INVALID);
-  assert(candidate_faces.size() > 0);
-  v_index = outer_vid;
-  A.resize(candidate_faces.size());
-  std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
+    assert(outer_vid != INVALID);
+    assert(candidate_faces.size() > 0);
+    v_index = outer_vid;
+    A.resize(candidate_faces.size());
+    std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
 }
 
 template<
@@ -85,94 +86,122 @@ template<
     typename DerivedA
     >
 IGL_INLINE void igl::outer_edge(
-       const Eigen::PlainObjectBase<DerivedV> & V,
-       const Eigen::PlainObjectBase<DerivedF> & F,
-       const Eigen::PlainObjectBase<DerivedI> & I,
-       IndexType & v1,
-       IndexType & v2,
-       Eigen::PlainObjectBase<DerivedA> & A) {
-   // Algorithm:
-   //    Find an outer vertex first.
-   //    Find the incident edge with largest slope when projected onto XY
-   //      plane.
-   //    If there is still a tie, break it using the projected slope onto ZX
-   //      plane.
-   //    If there is still a tie, then there are multiple overlapping edges,
-   //      which violates the precondition.
-   typedef typename DerivedV::Scalar Scalar;
-   typedef typename DerivedV::Index Index;
-   typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
-   typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
-   const size_t INVALID = std::numeric_limits<size_t>::max();
-
-   Index outer_vid;
-   Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
-   outer_vertex(V, F, I, outer_vid, candidate_faces);
-   const ScalarArray3& outer_v = V.row(outer_vid);
-   assert(candidate_faces.size() > 0);
-
-   auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
-   {
-     if (f[0] == vid) return 0;
-     if (f[1] == vid) return 1;
-     if (f[2] == vid) return 2;
-     assert(false);
-     return -1;
-   };
-
-   Scalar outer_slope_YX = 0;
-   Scalar outer_slope_ZX = 0;
-   size_t outer_opp_vid = INVALID;
-   bool infinite_slope_detected = false;
-   std::vector<Index> incident_faces;
-   auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
-     if (opp_vid == outer_opp_vid)
-     {
-       incident_faces.push_back(fid);
-       return;
-     }
-
-     const ScalarArray3 opp_v = V.row(opp_vid);
-     if (!infinite_slope_detected && outer_v[0] != opp_v[0]) 
-     {
-       // Finite slope
-       const ScalarArray3 diff = opp_v - outer_v;
-       const Scalar slope_YX = diff[1] / diff[0];
-       const Scalar slope_ZX = diff[2] / diff[0];
-       if (outer_opp_vid == INVALID ||
-               slope_YX > outer_slope_YX ||
-               (slope_YX == outer_slope_YX &&
-                slope_ZX > outer_slope_ZX)) {
-           outer_opp_vid = opp_vid;
-           outer_slope_YX = slope_YX;
-           outer_slope_ZX = slope_ZX;
-           incident_faces = {fid};
-       }
-     } else if (!infinite_slope_detected)
-     {
-       // Infinite slope
-       outer_opp_vid = opp_vid;
-       infinite_slope_detected = true;
-       incident_faces = {fid};
-     }
-   };
-
-   const size_t num_candidate_faces = candidate_faces.size();
-   for (size_t i=0; i<num_candidate_faces; i++) 
-   {
-     const Index fid = candidate_faces(i);
-     const IndexArray3& f = F.row(fid);
-     size_t id = get_vertex_index(f, outer_vid);
-     Index next_vid = f((id+1)%3);
-     Index prev_vid = f((id+2)%3);
-     check_and_update_outer_edge(next_vid, fid);
-     check_and_update_outer_edge(prev_vid, fid);
-   }
-
-   v1 = outer_vid;
-   v2 = outer_opp_vid;
-   A.resize(incident_faces.size());
-   std::copy(incident_faces.begin(), incident_faces.end(), A.data());
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & v1,
+        IndexType & v2,
+        Eigen::PlainObjectBase<DerivedA> & A) {
+    // Algorithm:
+    //    Find an outer vertex first.
+    //    Find the incident edge with largest abs slope when projected onto XY plane.
+    //    If there is a tie, check the signed slope and use the positive one.
+    //    If there is still a tie, break it using the projected slope onto ZX plane.
+    //    If there is still a tie, again check the signed slope and use the positive one.
+    //    If there is still a tie, then there are multiple overlapping edges,
+    //    which violates the precondition.
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
+    typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
+    const Index INVALID = std::numeric_limits<Index>::max();
+
+    Index outer_vid;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
+    outer_vertex(V, F, I, outer_vid, candidate_faces);
+    const ScalarArray3& outer_v = V.row(outer_vid);
+    assert(candidate_faces.size() > 0);
+
+    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
+    {
+        if (f[0] == vid) return 0;
+        if (f[1] == vid) return 1;
+        if (f[2] == vid) return 2;
+        assert(false);
+        return -1;
+    };
+
+    auto unsigned_value = [](Scalar v) -> Scalar {
+        if (v < 0) return v * -1;
+        else return v;
+    };
+
+    Scalar outer_slope_YX = 0;
+    Scalar outer_slope_ZX = 0;
+    Index outer_opp_vid = INVALID;
+    bool infinite_slope_detected = false;
+    std::vector<Index> incident_faces;
+    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
+        if (opp_vid == outer_opp_vid)
+        {
+            incident_faces.push_back(fid);
+            return;
+        }
+
+        const ScalarArray3 opp_v = V.row(opp_vid);
+        if (!infinite_slope_detected && outer_v[0] != opp_v[0])
+        {
+            // Finite slope
+            const ScalarArray3 diff = opp_v - outer_v;
+            const Scalar slope_YX = diff[1] / diff[0];
+            const Scalar slope_ZX = diff[2] / diff[0];
+            const Scalar u_slope_YX = unsigned_value(slope_YX);
+            const Scalar u_slope_ZX = unsigned_value(slope_ZX);
+            bool update = false;
+            if (outer_opp_vid == INVALID) {
+                update = true;
+            } else {
+                const Scalar u_outer_slope_YX = unsigned_value(outer_slope_YX);
+                if (u_slope_YX > u_outer_slope_YX) {
+                    update = true;
+                } else if (u_slope_YX == u_outer_slope_YX &&
+                        slope_YX > outer_slope_YX) {
+                    update = true;
+                } else if (slope_YX == outer_slope_YX) {
+                    const Scalar u_outer_slope_ZX =
+                        unsigned_value(outer_slope_ZX);
+                    if (u_slope_ZX > u_outer_slope_ZX) {
+                        update = true;
+                    } else if (u_slope_ZX == u_outer_slope_ZX &&
+                            slope_ZX > outer_slope_ZX) {
+                        update = true;
+                    } else if (slope_ZX == u_outer_slope_ZX) {
+                        assert(false);
+                    }
+                }
+            }
+
+            if (update) {
+                outer_opp_vid = opp_vid;
+                outer_slope_YX = slope_YX;
+                outer_slope_ZX = slope_ZX;
+                incident_faces = {fid};
+            }
+        } else if (!infinite_slope_detected)
+        {
+            // Infinite slope
+            outer_opp_vid = opp_vid;
+            infinite_slope_detected = true;
+            incident_faces = {fid};
+        }
+    };
+
+    const size_t num_candidate_faces = candidate_faces.size();
+    for (size_t i=0; i<num_candidate_faces; i++)
+    {
+        const Index fid = candidate_faces(i);
+        const IndexArray3& f = F.row(fid);
+        size_t id = get_vertex_index(f, outer_vid);
+        Index next_vid = f((id+1)%3);
+        Index prev_vid = f((id+2)%3);
+        check_and_update_outer_edge(next_vid, fid);
+        check_and_update_outer_edge(prev_vid, fid);
+    }
+
+    v1 = outer_vid;
+    v2 = outer_opp_vid;
+    A.resize(incident_faces.size());
+    std::copy(incident_faces.begin(), incident_faces.end(), A.data());
 }
 
 template<

include/igl/remove_unreferenced.cpp

@@ -45,7 +45,8 @@ IGL_INLINE void igl::remove_unreferenced(
   const size_t n = V.rows();
   remove_unreferenced(n,F,I,J);
   NF = F;
-  for_each(NF.data(),NF.data()+NF.size(),[&I](int & a){a=I(a);});
+  for_each(NF.data(),NF.data()+NF.size(),
+    [&I](typename DerivedNF::Scalar & a){a=I(a);});
   slice(V,J,1,NV);
 }
 

include/igl/slice.cpp

@@ -134,14 +134,14 @@ IGL_INLINE void igl::slice(
 #endif
 }
 
-template <typename MatX, typename MatY>
+template <typename MatX, typename DerivedR, typename MatY>
 IGL_INLINE void igl::slice(
   const MatX& X,
-  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+  const Eigen::PlainObjectBase<DerivedR> & R,
   const int dim,
   MatY& Y)
 {
-  Eigen::VectorXi C;
+  Eigen::Matrix<typename DerivedR::Scalar,Eigen::Dynamic,1> C;
   switch(dim)
   {
     case 1:
@@ -168,11 +168,15 @@ IGL_INLINE void igl::slice(
   }
 }
 
-template <typename DerivedX, typename DerivedY>
+template <
+  typename DerivedX, 
+  typename DerivedR,
+  typename DerivedC,
+  typename DerivedY>
 IGL_INLINE void igl::slice(
   const Eigen::PlainObjectBase<DerivedX> & X,
-  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
-  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+  const Eigen::PlainObjectBase<DerivedR> & R,
+  const Eigen::PlainObjectBase<DerivedC> & C,
   Eigen::PlainObjectBase<DerivedY> & Y)
 {
 #ifndef NDEBUG
@@ -244,20 +248,20 @@ IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::slice<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > igl::slice<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > igl::slice<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
-template void igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
+template void igl::slice<std::complex<double>, std::complex<double> >(Eigen::SparseMatrix<std::complex<double>, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<std::complex<double>, 0, int>&);
+template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template void igl::slice<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::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::slice<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::slice<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::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 template void igl::slice<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-template void igl::slice<std::complex<double>, std::complex<double> >(Eigen::SparseMatrix<std::complex<double>, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<std::complex<double>, 0, int>&);
-template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/slice.h

@@ -23,7 +23,9 @@ namespace igl
   //   Y  #R by #C matrix
   //
   // See also: slice_mask
-  template <typename TX, typename TY>
+  template <
+    typename TX, 
+    typename TY>
   IGL_INLINE void slice(
     const Eigen::SparseMatrix<TX>& X,
     const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
@@ -35,17 +37,24 @@ namespace igl
   //   dim  dimension to slice in 1 or 2, dim=1 --> X(R,:), dim=2 --> X(:,R)
   //
   // Note: For now this is just a cheap wrapper.
-  template <typename MatX, typename MatY>
+  template <
+    typename MatX, 
+    typename DerivedR,
+    typename MatY>
   IGL_INLINE void slice(
     const MatX& X,
-    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+    const Eigen::PlainObjectBase<DerivedR> & R,
     const int dim,
     MatY& Y);
-  template <typename DerivedX, typename DerivedY>
+  template <
+    typename DerivedX, 
+    typename DerivedR, 
+    typename DerivedC, 
+    typename DerivedY>
   IGL_INLINE void slice(
     const Eigen::PlainObjectBase<DerivedX> & X,
-    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
-    const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+    const Eigen::PlainObjectBase<DerivedR> & R,
+    const Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedY> & Y);
 
   template <typename DerivedX, typename DerivedY>

tutorial/509_Planarization/main.cpp

@@ -90,19 +90,19 @@ int main(int argc, char *argv[])
   FQCtri.resize(2*FQC.rows(), 3);
   FQCtri <<  FQC.col(0),FQC.col(1),FQC.col(2),
              FQC.col(2),FQC.col(3),FQC.col(0);
-  igl::slice( VQC, FQC.col(0), 1, PQC0);
-  igl::slice( VQC, FQC.col(1), 1, PQC1);
-  igl::slice( VQC, FQC.col(2), 1, PQC2);
-  igl::slice( VQC, FQC.col(3), 1, PQC3);
+  igl::slice( VQC, FQC.col(0).eval(), 1, PQC0);
+  igl::slice( VQC, FQC.col(1).eval(), 1, PQC1);
+  igl::slice( VQC, FQC.col(2).eval(), 1, PQC2);
+  igl::slice( VQC, FQC.col(3).eval(), 1, PQC3);
 
   // Planarize it
   igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCplan);
 
   // Convert the planarized mesh to triangles
-  igl::slice( VQCplan, FQC.col(0), 1, PQC0plan);
-  igl::slice( VQCplan, FQC.col(1), 1, PQC1plan);
-  igl::slice( VQCplan, FQC.col(2), 1, PQC2plan);
-  igl::slice( VQCplan, FQC.col(3), 1, PQC3plan);
+  igl::slice( VQCplan, FQC.col(0).eval(), 1, PQC0plan);
+  igl::slice( VQCplan, FQC.col(1).eval(), 1, PQC1plan);
+  igl::slice( VQCplan, FQC.col(2).eval(), 1, PQC2plan);
+  igl::slice( VQCplan, FQC.col(3).eval(), 1, PQC3plan);
 
   // Launch the viewer
   igl::viewer::Viewer viewer;

tutorial/609_Boolean/main.cpp

@@ -1,6 +1,6 @@
 #include <igl/readOFF.h>
 #define IGL_NO_CORK
-#undef IGL_STATIC_LIBRARY
+//#undef IGL_STATIC_LIBRARY
 #include <igl/boolean/mesh_boolean.h>
 #include <igl/viewer/Viewer.h>
 
@@ -8,7 +8,7 @@
 #include <iostream>
 
 Eigen::MatrixXd VA,VB,VC;
-Eigen::VectorXi J;
+Eigen::VectorXi J,I;
 Eigen::MatrixXi FA,FB,FC;
 igl::boolean::MeshBooleanType boolean_type(
   igl::boolean::MESH_BOOLEAN_TYPE_UNION);
@@ -24,7 +24,7 @@ const char * MESH_BOOLEAN_TYPE_NAMES[] =
 
 void update(igl::viewer::Viewer &viewer)
 {
-  igl::boolean::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC,J);
+  igl::boolean::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC,J,I);
   Eigen::MatrixXd C(FC.rows(),3);
   for(size_t f = 0;f<C.rows();f++)
   {