Merge branch 'master' of github.com:libigl/libigl into gh-pages

Former-commit-id: 03403d94f19182b4a1977966d6bd7bb3ba21b0a4

README.md

@@ -191,19 +191,22 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
  - [Adobe Research](http://www.adobe.com/technology/)  
+ - [Mesh](http://meshconsultants.ca/), consultants, Canada
  - [Pixar Research](http://graphics.pixar.com/research/)
  - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.
  - Columbia University, [Columbia Computer Graphics Group](http://www.cs.columbia.edu/cg/), USA
  - [Cornell University](http://www.graphics.cornell.edu/), USA
+ - [Czech Technical University in Prague](http://dcgi.felk.cvut.cz/), Czech
  - EPF Lausanne, [Computer Graphics and Geometry Laboratory](http://lgg.epfl.ch/people.php), Switzerland
  - ETH Zurich, [Interactive Geometry Lab](http://igl.ethz.ch/) and [Advanced Technologies Lab](http://ait.inf.ethz.ch/), Swizterland
  - George Mason University, [CraGL](http://cs.gmu.edu/~ygingold/), USA
- - [Hong Kong University of Science and Technology](http://www.ust.hk/), USA
+ - [Hong Kong University of Science and Technology](http://www.ust.hk/), Hong Kong
  - [National Institute of Informatics](http://www.nii.ac.jp/en/), Japan
  - New York University, [Media Research Lab](http://mrl.nyu.edu/), USA
  - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
  - [TU Berlin](https://www.cg.tu-berlin.de), Germany
  - [TU Delft](http://www.tudelft.nl/en/), Netherlands
+ - [TU Wien](https://www.tuwien.ac.at/en/tuwien_home/), Austria
  - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
  - [Universidade Federal de Santa Catarina](http://mtm.ufsc.br/~leo/), Brazil
  - [University College London](http://vecg.cs.ucl.ac.uk/), England
@@ -214,6 +217,7 @@ few labs/companies/institutions using libigl:
  - [University of Victoria](https://www.csc.uvic.ca/Research/graphics/), Canada
  - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
  - [Université Toulouse III Paul Sabatier](http://www.univ-tlse3.fr/), France
+ - [Zhejiang University](http://www.math.zju.edu.cn/cagd/), China
 
 
 ## Contact

include/igl/AABB.h

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

include/igl/ambient_occlusion.cpp

@@ -32,9 +32,9 @@ IGL_INLINE void igl::ambient_occlusion(
   S.resize(n,1);
   VectorXi hits = VectorXi::Zero(n,1);
   // Embree seems to be parallel when constructing but not when tracing rays
-#pragma omp parallel for
   const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
   // loop over mesh vertices
+#pragma omp parallel for
   for(int p = 0;p<n;p++)
   {
     const Vector3f origin = P.row(p).template cast<float>();

include/igl/barycentric_coordinates.cpp

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

include/igl/barycentric_coordinates.h

@@ -39,15 +39,13 @@ namespace igl
   // Compute barycentric coordinates in a triangle
   //
   // Inputs:
-  //   P  #P by 2 Query points in 2d
-  //   A  #P by 2 Triangle corners in 2d
-  //   B  #P by 2 Triangle corners in 2d
-  //   C  #P by 2 Triangle corners in 2d
+  //   P  #P by dim Query points
+  //   A  #P by dim Triangle corners
+  //   B  #P by dim Triangle corners
+  //   C  #P by dim Triangle corners
   // Outputs:
-  //   L  #P by e list of barycentric coordinates
+  //   L  #P by 3 list of barycentric coordinates
   //   
-  // Known bugs: this code is not tested (and probably will not work) for
-  // triangles and queries in 3D even if the query lives in/on the triangle.
   template <
     typename DerivedP,
     typename DerivedA,
@@ -55,10 +53,10 @@ namespace igl
     typename DerivedC,
     typename DerivedL>
   IGL_INLINE void barycentric_coordinates(
-    const Eigen::PlainObjectBase<DerivedP> & P,
-    const Eigen::PlainObjectBase<DerivedA> & A,
-    const Eigen::PlainObjectBase<DerivedB> & B,
-    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const Eigen::MatrixBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedL> & L);
 
 }

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -16,11 +16,13 @@
 #include "../../unique_simplices.h"
 #include "../../slice.h"
 #include "../../resolve_duplicated_faces.h"
+#include "../../get_seconds.h"
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <algorithm>
 
 //#define MESH_BOOLEAN_TIMING
+//#define DOUBLE_CHECK_EXACT_OUTPUT
 
 template <
   typename DerivedVA,
@@ -33,7 +35,7 @@ template <
   typename DerivedVC,
   typename DerivedFC,
   typename DerivedJ>
-IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
+IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
     const Eigen::PlainObjectBase<DerivedVA> & VA,
     const Eigen::PlainObjectBase<DerivedFA> & FA,
     const Eigen::PlainObjectBase<DerivedVB> & VB,
@@ -43,7 +45,8 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     const ResolveFunc& resolve_fun,
     Eigen::PlainObjectBase<DerivedVC > & VC,
     Eigen::PlainObjectBase<DerivedFC > & FC,
-    Eigen::PlainObjectBase<DerivedJ > & J) {
+    Eigen::PlainObjectBase<DerivedJ > & J) 
+{
 
 #ifdef MESH_BOOLEAN_TIMING
   const auto & tictoc = []() -> double
@@ -82,23 +85,6 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
       DerivedFC FF(FA.rows() + FB.rows(), 3);
       VV << VA, VB;
       FF << FA, FB.array() + VA.rows();
-      //// Handle annoying empty cases
-      //if(VA.size()>0)
-      //{
-      //  VV<<VA;
-      //}
-      //if(VB.size()>0)
-      //{
-      //  VV<<VB;
-      //}
-      //if(FA.size()>0)
-      //{
-      //  FF<<FA;
-      //}
-      //if(FB.size()>0)
-      //{
-      //  FF<<FB.array()+VA.rows();
-      //}
       resolve_fun(VV, FF, V, F, CJ);
   }
 #ifdef MESH_BOOLEAN_TIMING
@@ -111,13 +97,18 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
   Eigen::VectorXi labels(num_faces);
   std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
       [&](int i) { return i<FA.rows() ? 0:1; });
-  if (num_faces > 0) {
-    igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
-  } else {
+  bool valid = true;
+  if (num_faces > 0) 
+  {
+    valid = valid & 
+      igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
+  } else 
+  {
     W.resize(0, 4);
   }
   assert((size_t)W.rows() == num_faces);
-  if (W.cols() == 2) {
+  if (W.cols() == 2) 
+  {
     assert(FB.rows() == 0);
     Eigen::MatrixXi W_tmp(num_faces, 4);
     W_tmp << W, Eigen::MatrixXi::Zero(num_faces, 2);
@@ -131,7 +122,8 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
 
   // Compute resulting winding number.
   Eigen::MatrixXi Wr(num_faces, 2);
-  for (size_t i=0; i<num_faces; i++) {
+  for (size_t i=0; i<num_faces; i++) 
+  {
     Eigen::MatrixXi w_out(1,2), w_in(1,2);
     w_out << W(i,0), W(i,2);
     w_in  << W(i,1), W(i,3);
@@ -143,18 +135,22 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
 #endif
 
   // Extract boundary separating inside from outside.
-  auto index_to_signed_index = [&](size_t i, bool ori) -> int{
+  auto index_to_signed_index = [&](size_t i, bool ori) -> int
+  {
     return (i+1)*(ori?1:-1);
   };
   //auto signed_index_to_index = [&](int i) -> size_t {
   //    return abs(i) - 1;
   //};
   std::vector<int> selected;
-  for(size_t i=0; i<num_faces; i++) {
+  for(size_t i=0; i<num_faces; i++) 
+  {
     auto should_keep = keep(Wr(i,0), Wr(i,1));
-    if (should_keep > 0) {
+    if (should_keep > 0) 
+    {
       selected.push_back(index_to_signed_index(i, true));
-    } else if (should_keep < 0) {
+    } else if (should_keep < 0) 
+    {
       selected.push_back(index_to_signed_index(i, false));
     }
   }
@@ -162,11 +158,14 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
   const size_t num_selected = selected.size();
   DerivedFC kept_faces(num_selected, 3);
   DerivedJ  kept_face_indices(num_selected, 1);
-  for (size_t i=0; i<num_selected; i++) {
+  for (size_t i=0; i<num_selected; i++) 
+  {
     size_t idx = abs(selected[i]) - 1;
-    if (selected[i] > 0) {
+    if (selected[i] > 0) 
+    {
       kept_faces.row(i) = F.row(idx);
-    } else {
+    } else 
+    {
       kept_faces.row(i) = F.row(idx).reverse();
     }
     kept_face_indices(i, 0) = CJ[idx];
@@ -182,6 +181,31 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     igl::resolve_duplicated_faces(kept_faces, G, JJ);
     igl::slice(kept_face_indices, JJ, 1, J);
 
+#ifdef DOUBLE_CHECK_EXACT_OUTPUT
+    {
+      // Sanity check on exact output.
+      igl::copyleft::cgal::RemeshSelfIntersectionsParam params;
+      params.detect_only = true;
+      params.first_only = true;
+      MatrixXES dummy_VV;
+      DerivedFC dummy_FF, dummy_IF;
+      Eigen::VectorXi dummy_J, dummy_IM;
+      igl::copyleft::cgal::SelfIntersectMesh<
+        Kernel,
+        MatrixXES, DerivedFC,
+        MatrixXES, DerivedFC,
+        DerivedFC,
+        Eigen::VectorXi,
+        Eigen::VectorXi
+      > checker(V, G, params,
+          dummy_VV, dummy_FF, dummy_IF, dummy_J, dummy_IM);
+      if (checker.count != 0) 
+      {
+        throw "Self-intersection not fully resolved.";
+      }
+    }
+#endif
+
     MatrixX3S Vs(V.rows(), V.cols());
     for (size_t i=0; i<(size_t)V.rows(); i++)
     {
@@ -196,6 +220,7 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
 #ifdef MESH_BOOLEAN_TIMING
   log_time("clean_up");
 #endif
+  return valid;
 }
 
 template <
@@ -207,7 +232,7 @@ template <
   typename DerivedVC,
   typename DerivedFC,
   typename DerivedJ>
-IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
+IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
     const Eigen::PlainObjectBase<DerivedVA > & VA,
     const Eigen::PlainObjectBase<DerivedFA > & FA,
     const Eigen::PlainObjectBase<DerivedVB > & VB,
@@ -218,43 +243,32 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     Eigen::PlainObjectBase<DerivedFC > & FC,
     Eigen::PlainObjectBase<DerivedJ > & J)
 {
-  //typedef CGAL::Epeck Kernel;
-  //typedef Kernel::FT ExactScalar;
-  //typedef Eigen::Matrix<
-  //  ExactScalar,
-  //  Eigen::Dynamic,
-  //  Eigen::Dynamic,
-  //  DerivedVC::IsRowMajor> MatrixXES;
-
-  switch (type) {
+  switch (type) 
+  {
     case MESH_BOOLEAN_TYPE_UNION:
-      igl::copyleft::boolean::mesh_boolean(
+      return igl::copyleft::boolean::mesh_boolean(
           VA, FA, VB, FB, igl::copyleft::boolean::BinaryUnion(),
           igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-      break;
     case MESH_BOOLEAN_TYPE_INTERSECT:
-      igl::copyleft::boolean::mesh_boolean(
+      return igl::copyleft::boolean::mesh_boolean(
           VA, FA, VB, FB, igl::copyleft::boolean::BinaryIntersect(),
           igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-      break;
     case MESH_BOOLEAN_TYPE_MINUS:
-      igl::copyleft::boolean::mesh_boolean(
+      return igl::copyleft::boolean::mesh_boolean(
           VA, FA, VB, FB, igl::copyleft::boolean::BinaryMinus(),
           igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-      break;
     case MESH_BOOLEAN_TYPE_XOR:
-      igl::copyleft::boolean::mesh_boolean(
+      return igl::copyleft::boolean::mesh_boolean(
           VA, FA, VB, FB, igl::copyleft::boolean::BinaryXor(),
           igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-      break;
     case MESH_BOOLEAN_TYPE_RESOLVE:
       //op = binary_resolve();
-      igl::copyleft::boolean::mesh_boolean(
+      return igl::copyleft::boolean::mesh_boolean(
           VA, FA, VB, FB, igl::copyleft::boolean::BinaryResolve(),
           igl::copyleft::boolean::KeepAll(), resolve_func, VC, FC, J);
-      break;
     default:
-      throw std::runtime_error("Unsupported boolean type.");
+      assert(false && "Unsupported boolean type.");
+      return false;
   }
 }
 
@@ -266,7 +280,7 @@ template <
   typename DerivedVC,
   typename DerivedFC,
   typename DerivedJ>
-IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
+IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
   const Eigen::PlainObjectBase<DerivedVA > & VA,
   const Eigen::PlainObjectBase<DerivedFA > & FA,
   const Eigen::PlainObjectBase<DerivedVB > & VB,
@@ -318,29 +332,30 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
 }
 
 template <
-typename DerivedVA,
-typename DerivedFA,
-typename DerivedVB,
-typename DerivedFB,
-typename DerivedVC,
-typename DerivedFC>
-IGL_INLINE void igl::copyleft::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) {
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC>
+IGL_INLINE bool igl::copyleft::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;
   return igl::copyleft::boolean::mesh_boolean(VA,FA,VB,FB,type,VC,FC,J);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::copyleft::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::copyleft::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::copyleft::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::copyleft::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> >&);
-template void igl::copyleft::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::copyleft::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 bool igl::copyleft::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::copyleft::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::copyleft::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::copyleft::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> >&);
+template bool igl::copyleft::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::copyleft::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> >&);
 #undef IGL_STATIC_LIBRARY
 #include "../../remove_unreferenced.cpp"
 template void igl::remove_unreferenced<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<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/copyleft/boolean/mesh_boolean.h

@@ -39,6 +39,8 @@ namespace igl
       //    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
+      //  Returns true iff inputs induce a piecewise constant winding number
+      //    field
       //
       //  See also: mesh_boolean_cork, intersect_other,
       //  remesh_self_intersections
@@ -53,7 +55,7 @@ namespace igl
         typename DerivedVC,
         typename DerivedFC,
         typename DerivedJ>
-      IGL_INLINE void mesh_boolean(
+      IGL_INLINE bool mesh_boolean(
           const Eigen::PlainObjectBase<DerivedVA> & VA,
           const Eigen::PlainObjectBase<DerivedFA> & FA,
           const Eigen::PlainObjectBase<DerivedVB> & VB,
@@ -77,6 +79,8 @@ namespace igl
       //    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
+      //  Returns true if inputs induce a piecewise constant winding number
+      //    field and type is valid.
       //
       //  See also: mesh_boolean_cork, intersect_other,
       //  remesh_self_intersections
@@ -89,7 +93,7 @@ namespace igl
         typename DerivedVC,
         typename DerivedFC,
         typename DerivedJ>
-      IGL_INLINE void mesh_boolean(
+      IGL_INLINE bool mesh_boolean(
           const Eigen::PlainObjectBase<DerivedVA > & VA,
           const Eigen::PlainObjectBase<DerivedFA > & FA,
           const Eigen::PlainObjectBase<DerivedVB > & VB,
@@ -110,6 +114,8 @@ namespace igl
       //    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
+      //  Returns true if inputs induce a piecewise constant winding number
+      //  field and type is valid
       //
       //  See also: mesh_boolean_cork, intersect_other,
       //  remesh_self_intersections
@@ -121,7 +127,7 @@ namespace igl
         typename DerivedVC,
         typename DerivedFC,
         typename DerivedJ>
-      IGL_INLINE void mesh_boolean(
+      IGL_INLINE bool mesh_boolean(
         const Eigen::PlainObjectBase<DerivedVA > & VA,
         const Eigen::PlainObjectBase<DerivedFA > & FA,
         const Eigen::PlainObjectBase<DerivedVB > & VB,
@@ -140,6 +146,8 @@ namespace igl
       //  Outputs:
       //    VC  #VC by 3 list of vertex positions of boolean result mesh
       //    FC  #FC by 3 list of triangle indices into VC
+      //  Returns true ff inputs induce a piecewise constant winding number
+      //    field and type is valid
       template <
         typename DerivedVA,
         typename DerivedFA,
@@ -147,7 +155,7 @@ namespace igl
         typename DerivedFB,
         typename DerivedVC,
         typename DerivedFC>
-      IGL_INLINE void mesh_boolean(
+      IGL_INLINE bool mesh_boolean(
           const Eigen::PlainObjectBase<DerivedVA > & VA,
           const Eigen::PlainObjectBase<DerivedFA > & FA,
           const Eigen::PlainObjectBase<DerivedVB > & VB,

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -10,6 +10,7 @@
 
 #include "CGAL_includes.hpp"
 #include "RemeshSelfIntersectionsParam.h"
+#include "../../unique.h"
 
 #include <Eigen/Dense>
 #include <list>
@@ -423,7 +424,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     return;
   }
 
-  remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
+  remesh_intersections(V,F,T,offending,edge2faces,true,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");
@@ -702,8 +703,6 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
 {
   using namespace std;
 
-  {
-    std::lock_guard<std::mutex> guard(m_offending_lock);
   // must be co-planar
   if(
     A.supporting_plane() != B.supporting_plane() &&
@@ -711,7 +710,6 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     return false;
   }
-  }
   // Since A and B are non-degenerate the intersection must be a polygon
   // (triangle). Either
   //   - the vertex of A (B) opposite the shared edge of lies on B (A), or
@@ -865,22 +863,46 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   DerivedJ,
   DerivedIM>::process_intersecting_boxes()
 {
-  CGAL::Identity_transformation I;
   std::vector<std::mutex> triangle_locks(T.size());
-  auto process_chunk = [&](size_t first, size_t last) -> void{
+  std::vector<std::mutex> vertex_locks(V.rows());
+  std::mutex index_lock;
+  auto process_chunk = [&](const size_t first, const size_t last) -> void{
     assert(last >= first);
 
     for (size_t i=first; i<last; i++) {
-      const auto& box_pair = candidate_box_pairs[i];
-      const auto& a = box_pair.first;
-      const auto& b = box_pair.second;
-      Index fa = a.handle()-T.begin();
-      Index fb = b.handle()-T.begin();
+      Index fa=T.size(), fb=T.size();
+      {
+        // Before knowing which triangles are involved, we need to lock
+        // everything to prevent race condition in updating reference counters.
+        std::lock_guard<std::mutex> guard(index_lock);
+        const auto& box_pair = candidate_box_pairs[i];
+        const auto& a = box_pair.first;
+        const auto& b = box_pair.second;
+        fa = a.handle()-T.begin();
+        fb = b.handle()-T.begin();
+      }
+      assert(fa < T.size());
+      assert(fb < T.size());
 
+      // Lock triangles
       std::lock_guard<std::mutex> guard_A(triangle_locks[fa]);
       std::lock_guard<std::mutex> guard_B(triangle_locks[fb]);
-      const Triangle_3& A = *a.handle();
-      const Triangle_3& B = *b.handle();
+
+      // Lock vertices
+      std::list<std::lock_guard<std::mutex> > guard_vertices;
+      {
+        std::vector<typename DerivedF::Scalar> unique_vertices;
+        std::vector<size_t> tmp1, tmp2;
+        igl::unique({F(fa,0), F(fa,1), F(fa,2), F(fb,0), F(fb,1), F(fb,2)},
+            unique_vertices, tmp1, tmp2);
+        std::for_each(unique_vertices.begin(), unique_vertices.end(),
+            [&](const typename DerivedF::Scalar& vi) {
+            guard_vertices.emplace_back(vertex_locks[vi]);
+            });
+      }
+
+      const Triangle_3& A = T[fa];
+      const Triangle_3& B = T[fb];
 
       // Number of combinatorially shared vertices
       Index comb_shared_vertices = 0;
@@ -941,7 +963,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
         single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
       }else
       {
-        intersect(*a.handle(),*b.handle(),fa,fb);
+        intersect(A,B,fa,fb);
       }
     }
   };

include/igl/copyleft/cgal/closest_facet.cpp

@@ -12,13 +12,8 @@
 #include <stdexcept>
 #include <unordered_map>
 
-#include <CGAL/AABB_tree.h>
-#include <CGAL/AABB_traits.h>
-#include <CGAL/AABB_triangle_primitive.h>
-#include <CGAL/intersections.h>
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-
 #include "order_facets_around_edge.h"
+#include "submesh_aabb_tree.h"
 #include "../../vertex_triangle_adjacency.h"
 #include "../../writePLY.h"
 
@@ -39,7 +34,9 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     const std::vector<std::vector<uE2EType> >& uE2E,
     const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
     Eigen::PlainObjectBase<DerivedR>& R,
-    Eigen::PlainObjectBase<DerivedS>& S) {
+    Eigen::PlainObjectBase<DerivedS>& S)
+{
+
   typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
   typedef Kernel::Point_3 Point_3;
   typedef Kernel::Plane_3 Plane_3;
@@ -55,29 +52,64 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
         "Closest facet cannot be computed on empty mesh.");
   }
 
-  std::vector<std::vector<size_t> > VF;
-  std::vector<std::vector<size_t> > VFi;
+  std::vector<std::vector<size_t> > VF, VFi;
   igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
+  std::vector<bool> in_I;
+  std::vector<Triangle> triangles;
+  Tree tree;
+  submesh_aabb_tree(V,F,I,tree,triangles,in_I);
+
+  return closest_facet(
+    V,F,I,P,uE2E,EMAP,VF,VFi,tree,triangles,in_I,R,S);
+}
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename Kernel,
+  typename DerivedR,
+  typename DerivedS >
+IGL_INLINE void igl::copyleft::cgal::closest_facet(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedI>& I,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    const std::vector<std::vector<size_t> > & VF,
+    const std::vector<std::vector<size_t> > & VFi,
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<
+        Kernel, 
+        CGAL::AABB_triangle_primitive<
+          Kernel, typename std::vector<
+            typename Kernel::Triangle_3 >::iterator > > > & tree,
+    const std::vector<typename Kernel::Triangle_3 > & triangles,
+    const std::vector<bool> & in_I,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S)
+{
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename Kernel::Plane_3 Plane_3;
+  typedef typename Kernel::Segment_3 Segment_3;
+  typedef typename Kernel::Triangle_3 Triangle;
+  typedef typename std::vector<Triangle>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
 
-  std::vector<bool> in_I(F.rows(), false);
   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));
-    in_I[I(i,0)] = true;
-    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)));
-    if (triangles.back().is_degenerate()) {
-      throw std::runtime_error(
-          "Input facet components contains degenerated triangles");
-    }
+  if (F.rows() <= 0 || I.rows() <= 0) {
+    throw std::runtime_error(
+        "Closest facet cannot be computed on empty mesh.");
   }
-  Tree tree(triangles.begin(), triangles.end());
-  tree.accelerate_distance_queries();
 
-  auto on_the_positive_side = [&](size_t fid, const Point_3& p) -> bool{
+  auto on_the_positive_side = [&](size_t fid, const Point_3& p) -> bool
+  {
     const auto& f = F.row(fid).eval();
     Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
     Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
@@ -99,7 +131,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     return false;
   };
 
-  auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
+  auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool 
+  {
     const auto& f = F.row(fid);
     if      ((size_t)f[0] == s && (size_t)f[1] == d) return false;
     else if ((size_t)f[1] == s && (size_t)f[2] == d) return false;
@@ -143,23 +176,28 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
       size_t query_idx,
       const size_t s, const size_t d,
       size_t preferred_facet,
-      bool& orientation) -> size_t {
+      bool& orientation) -> size_t 
+  {
     Point_3 query_point(
-        P(query_idx, 0),
-        P(query_idx, 1),
-        P(query_idx, 2));
+      P(query_idx, 0),
+      P(query_idx, 1),
+      P(query_idx, 2));
 
     size_t corner_idx = std::numeric_limits<size_t>::max();
     if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 1)) ||
-        (s == F(preferred_facet, 1) && d == F(preferred_facet, 0))) {
+        (s == F(preferred_facet, 1) && d == F(preferred_facet, 0))) 
+    {
       corner_idx = 2;
     } else if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 2)) ||
-        (s == F(preferred_facet, 2) && d == F(preferred_facet, 0))) {
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 0))) 
+    {
       corner_idx = 1;
     } else if ((s == F(preferred_facet, 1) && d == F(preferred_facet, 2)) ||
-        (s == F(preferred_facet, 2) && d == F(preferred_facet, 1))) {
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 1))) 
+    {
       corner_idx = 0;
-    } else {
+    } else 
+    {
       std::cerr << "s: " << s << "\t d:" << d << std::endl;
       std::cerr << F.row(preferred_facet) << std::endl;
       throw std::runtime_error(
@@ -169,9 +207,11 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     auto ueid = EMAP(preferred_facet + corner_idx * F.rows());
     auto eids = uE2E[ueid];
     std::vector<size_t> intersected_face_indices;
-    for (auto eid : eids) {
+    for (auto eid : eids) 
+    {
       const size_t fid = eid % F.rows();
-      if (in_I[fid]) {
+      if (in_I[fid]) 
+      {
         intersected_face_indices.push_back(fid);
       }
     }
@@ -188,7 +228,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
         });
 
     assert(num_intersected_faces >= 1);
-    if (num_intersected_faces == 1) {
+    if (num_intersected_faces == 1) 
+    {
       // The edge must be a boundary edge.  Thus, the orientation can be
       // simply determined by checking if the query point is on the
       // positive side of the facet.
@@ -200,8 +241,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     Eigen::VectorXi order;
     DerivedP pivot = P.row(query_idx).eval();
     igl::copyleft::cgal::order_facets_around_edge(V, F, s, d,
-        intersected_face_signed_indices,
-        pivot, order);
+      intersected_face_signed_indices,
+      pivot, order);
 
     // Although first and last are equivalent, make the choice based on
     // preferred_facet.
@@ -381,7 +422,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     bool fid_ori = false;
 
     // Gether all facets sharing the closest point.
-    std::vector<Tree::Primitive_id> intersected_faces;
+    typename std::vector<typename Tree::Primitive_id> intersected_faces;
     tree.all_intersected_primitives(Segment_3(closest_point, query),
         std::back_inserter(intersected_faces));
     const size_t num_intersected_faces = intersected_faces.size();
@@ -389,7 +430,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     std::transform(intersected_faces.begin(),
         intersected_faces.end(),
         intersected_face_indices.begin(),
-        [&](const Tree::Primitive_id& itr) -> size_t
+        [&](const typename Tree::Primitive_id& itr) -> size_t
         { return I(itr-triangles.begin(), 0); });
 
     size_t element_index;

include/igl/copyleft/cgal/closest_facet.h

@@ -13,6 +13,12 @@
 #include <Eigen/Core>
 #include <vector>
 
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 namespace igl
 {
   namespace copyleft
@@ -32,14 +38,14 @@ namespace igl
       //   S  #P list of bools indicating on which side of the closest facet
       //      each query point lies.
       template<
-          typename DerivedV,
-          typename DerivedF,
-          typename DerivedI,
-          typename DerivedP,
-          typename uE2EType,
-          typename DerivedEMAP,
-          typename DerivedR,
-          typename DerivedS >
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename DerivedP,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename DerivedR,
+        typename DerivedS >
       IGL_INLINE void closest_facet(
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
@@ -49,23 +55,51 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
-
       template<
-          typename DerivedV,
-          typename DerivedF,
-          typename DerivedP,
-          typename uE2EType,
-          typename DerivedEMAP,
-          typename DerivedR,
-          typename DerivedS >
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename DerivedR,
+        typename DerivedS >
       IGL_INLINE void closest_facet(
-              const Eigen::PlainObjectBase<DerivedV>& V,
-              const Eigen::PlainObjectBase<DerivedF>& F,
-              const Eigen::PlainObjectBase<DerivedP>& P,
-              const std::vector<std::vector<uE2EType> >& uE2E,
-              const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-              Eigen::PlainObjectBase<DerivedR>& R,
-              Eigen::PlainObjectBase<DerivedS>& S);
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        Eigen::PlainObjectBase<DerivedR>& R,
+        Eigen::PlainObjectBase<DerivedS>& S);
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename DerivedP,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename Kernel,
+        typename DerivedR,
+        typename DerivedS >
+      IGL_INLINE void closest_facet(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const Eigen::PlainObjectBase<DerivedI>& I,
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          const std::vector<std::vector<uE2EType> >& uE2E,
+          const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+          const std::vector<std::vector<size_t> > & VF,
+          const std::vector<std::vector<size_t> > & VFi,
+          const CGAL::AABB_tree<
+            CGAL::AABB_traits<
+              Kernel, 
+              CGAL::AABB_triangle_primitive<
+                Kernel, typename std::vector<
+                  typename Kernel::Triangle_3 >::iterator > > > & tree,
+          const std::vector<typename Kernel::Triangle_3 > & triangles,
+          const std::vector<bool> & in_I,
+          Eigen::PlainObjectBase<DerivedR>& R,
+          Eigen::PlainObjectBase<DerivedS>& S);
     }
   }
 }

include/igl/copyleft/cgal/extract_cells.cpp

@@ -7,17 +7,28 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 //
 #include "extract_cells.h"
+#include "closest_facet.h"
+#include "order_facets_around_edge.h"
+#include "outer_facet.h"
+#include "submesh_aabb_tree.h"
 #include "../../extract_manifold_patches.h"
 #include "../../facet_components.h"
+#include "../../get_seconds.h"
 #include "../../triangle_triangle_adjacency.h"
 #include "../../unique_edge_map.h"
-#include "../../get_seconds.h"
-#include "closest_facet.h"
-#include "order_facets_around_edge.h"
-#include "outer_facet.h"
+#include "../../vertex_triangle_adjacency.h"
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
+#include <iostream>
 #include <vector>
 #include <queue>
+#include <map>
+#include <set>
 
 //#define EXTRACT_CELLS_DEBUG
 
@@ -72,6 +83,17 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
   const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
   Eigen::PlainObjectBase<DerivedC>& cells) 
 {
+
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+  typedef Kernel::Point_3 Point_3;
+  typedef Kernel::Plane_3 Plane_3;
+  typedef Kernel::Segment_3 Segment_3;
+  typedef Kernel::Triangle_3 Triangle;
+  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;
+
 #ifdef EXTRACT_CELLS_DEBUG
   const auto & tictoc = []() -> double
   {
@@ -118,11 +140,25 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     components[C[i]].push_back(i);
   }
   // Convert vector lists to Eigen lists...
+  // and precompute data-structures for each component
+  std::vector<std::vector<size_t> > VF,VFi;
+  igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
   std::vector<Eigen::VectorXi> Is(num_components);
+  std::vector<
+    CGAL::AABB_tree<
+      CGAL::AABB_traits<
+        Kernel, 
+        CGAL::AABB_triangle_primitive<
+          Kernel, std::vector<
+            Kernel::Triangle_3 >::iterator > > > > trees(num_components);
+  std::vector< std::vector<Kernel::Triangle_3 > > 
+    triangle_lists(num_components);
+  std::vector<std::vector<bool> > in_Is(num_components);
   for (size_t i=0; i<num_components; i++)
   {
     Is[i].resize(components[i].size());
     std::copy(components[i].begin(), components[i].end(),Is[i].data());
+    submesh_aabb_tree(V,F,Is[i],trees[i],triangle_lists[i],in_Is[i]);
   }
 
   // Find outer facets, their orientations and cells for each component
@@ -216,9 +252,25 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
       // Gather closest facets in ith component to each query point and their
       // orientations
       const auto& I = Is[i];
+      const auto& tree = trees[i];
+      const auto& in_I = in_Is[i];
+      const auto& triangles = triangle_lists[i];
+
       Eigen::VectorXi closest_facets, closest_facet_orientations;
-      closest_facet(V, F, I, queries,
-        uE2E, EMAP, closest_facets, closest_facet_orientations);
+      closest_facet(
+        V,
+        F, 
+        I, 
+        queries,
+        uE2E, 
+        EMAP, 
+        VF,
+        VFi,
+        tree,
+        triangles,
+        in_I,
+        closest_facets, 
+        closest_facet_orientations);
       // Loop over all candidates
       for (size_t j=0; j<num_candidate_comps; j++)
       {
@@ -366,186 +418,107 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells_single_component(
   const size_t num_unique_edges = uE.rows();
   const size_t num_patches = P.maxCoeff() + 1;
 
-  // patch_edge_adj[p] --> list {e,f,g,...} such that p is a patch id and
-  //   e,f,g etc. are edge indices into 
-  std::vector<std::vector<size_t> > patch_edge_adj(num_patches);
-  // orders[u] --> J  where u is an index into unique edges uE and J is a
-  //   #adjacent-faces list of face-edge indices into F*3 sorted cyclicaly around
-  //   edge u.
-  std::vector<Eigen::VectorXi> orders(num_unique_edges);
-  // orientations[u] ---> list {f1,f2, ...} where u is an index into unique edges uE
-  //   and points to #adj-facets long list of flags whether faces are oriented
-  //   to point their normals clockwise around edge when looking along the
-  //   edge.
-  std::vector<std::vector<bool> > orientations(num_unique_edges);
-  // Loop over unique edges
-  for (size_t i=0; i<num_unique_edges; i++) 
-  {
+  // Build patch-patch adjacency list.
+  std::vector<std::map<size_t, size_t> > patch_adj(num_patches);
+  for (size_t i=0; i<num_unique_edges; i++) {
     const size_t s = uE(i,0);
     const size_t d = uE(i,1);
     const auto adj_faces = uE2E[i];
-    // If non-manifold (more than two incident faces)
-    if (adj_faces.size() > 2) 
-    {
-      // signed_adj_faces[a] --> sid  where a is an index into adj_faces
-      // (list of face edges on {s,d}) and sid is a signed index for resolve
-      // co-planar duplicates consistently (i.e. using simulation of
-      // simplicity).
+    const size_t num_adj_faces = adj_faces.size();
+    if (num_adj_faces > 2) {
+      for (size_t j=0; j<num_adj_faces; j++) {
+        const size_t patch_j = P[edge_index_to_face_index(adj_faces[j])];
+        for (size_t k=j+1; k<num_adj_faces; k++) {
+          const size_t patch_k = P[edge_index_to_face_index(adj_faces[k])];
+          if (patch_adj[patch_j].find(patch_k) == patch_adj[patch_j].end()) {
+            patch_adj[patch_j].insert({patch_k, i});
+          }
+          if (patch_adj[patch_k].find(patch_j) == patch_adj[patch_k].end()) {
+            patch_adj[patch_k].insert({patch_j, i});
+          }
+        }
+      }
+    }
+  }
+
+
+  const int INVALID = std::numeric_limits<int>::max();
+  std::vector<size_t> cell_labels(num_patches * 2);
+  for (size_t i=0; i<num_patches; i++) cell_labels[i] = i;
+  std::vector<std::set<size_t> > equivalent_cells(num_patches*2);
+  std::vector<bool> processed(num_unique_edges, false);
+
+  size_t label_count=0;
+  for (size_t i=0; i<num_patches; i++) {
+    for (const auto& entry : patch_adj[i]) {
+      const size_t neighbor_patch = entry.first;
+      const size_t uei = entry.second;
+      if (processed[uei]) continue;
+      processed[uei] = true;
+
+      const auto& adj_faces = uE2E[uei];
+      const size_t num_adj_faces = adj_faces.size();
+      assert(num_adj_faces > 2);
+
+      const size_t s = uE(uei,0);
+      const size_t d = uE(uei,1);
+
       std::vector<int> signed_adj_faces;
-      for (auto ei : adj_faces)
+      for (auto ej : adj_faces)
       {
-        const size_t fid = edge_index_to_face_index(ei);
+        const size_t fid = edge_index_to_face_index(ej);
         bool cons = is_consistent(fid, s, d);
         signed_adj_faces.push_back((fid+1)*(cons ? 1:-1));
       }
       {
         // Sort adjacent faces cyclically around {s,d}
-        auto& order = orders[i];
+        Eigen::VectorXi order;
         // order[f] will reveal the order of face f in signed_adj_faces
         order_facets_around_edge(V, F, s, d, signed_adj_faces, order);
-        // Determine if each facet is oriented to point its normal clockwise or
-        // not around the {s,d} (normals are not explicitly computed/used)
-        auto& orientation = orientations[i];
-        orientation.resize(order.size());
-        std::transform(
-          order.data(), 
-          order.data() + order.size(),
-          orientation.begin(), 
-          [&signed_adj_faces](const int i){ return signed_adj_faces[i] > 0;});
-        // re-index order from adjacent faces to full face list. Now order
-        // indexes F directly
-        std::transform(
-          order.data(), 
-          order.data() + order.size(),
-          order.data(), 
-          [&adj_faces](const int index){ return adj_faces[index];});
-      }
-      // loop over adjacent faces, remember that patch is adjacent to this edge
-      for(const auto & ei : adj_faces)
-      {
-        const size_t fid = edge_index_to_face_index(ei);
-        patch_edge_adj[P[fid]].push_back(ei);
+        for (size_t j=0; j<num_adj_faces; j++) {
+          const size_t curr_idx = j;
+          const size_t next_idx = (j+1)%num_adj_faces;
+          const size_t curr_patch_idx =
+            P[edge_index_to_face_index(adj_faces[order[curr_idx]])];
+          const size_t next_patch_idx =
+            P[edge_index_to_face_index(adj_faces[order[next_idx]])];
+          const bool curr_cons = signed_adj_faces[order[curr_idx]] > 0;
+          const bool next_cons = signed_adj_faces[order[next_idx]] > 0;
+          const size_t curr_cell_idx = curr_patch_idx*2 + (curr_cons?0:1);
+          const size_t next_cell_idx = next_patch_idx*2 + (next_cons?1:0);
+          equivalent_cells[curr_cell_idx].insert(next_cell_idx);
+          equivalent_cells[next_cell_idx].insert(curr_cell_idx);
+        }
       }
     }
   }
 
-  // Initialize all patch-to-cell indices as "invalid" (unlabeled)
-  const int INVALID = std::numeric_limits<int>::max();
+  size_t count=0;
   cells.resize(num_patches, 2);
   cells.setConstant(INVALID);
-
-  // Given a "seed" patch id, a cell id, and which side of the patch that cell
-  // lies, identify all other patches bounding this cell (and tell them that
-  // they do)
-  //
-  // Inputs:
-  //   seed_patch_id  index into patches
-  //   cell_idx  index into cells
-  //   seed_patch_side   0 or 1 depending on whether cell_idx is on left or
-  //     right side of seed_patch_id 
-  //   cells  #cells by 2 list of current assignment of cells incident on each
-  //     side of a patch
-  // Outputs:
-  //   cells  udpated (see input)
-  // 
-  const auto & peel_cell_bd = 
-    [&P,&patch_edge_adj,&EMAP,&orders,&orientations,&num_faces](
-    const size_t seed_patch_id, 
-    const short seed_patch_side, 
-    const size_t cell_idx,
-    Eigen::PlainObjectBase<DerivedC>& cells)
-  {
-    typedef std::pair<size_t, short> PatchSide;
-    // Initialize a queue of {p,side} patch id and patch side pairs to BFS over
-    // all patches
-    std::queue<PatchSide> Q;
-    Q.emplace(seed_patch_id, seed_patch_side);
-    // assign cell id of seed patch on appropriate side
-    cells(seed_patch_id, seed_patch_side) = cell_idx;
-    while (!Q.empty())
-    {
-      // Pop patch from Q
-      const auto entry = Q.front();
+  const auto extract_equivalent_cells = [&](size_t i) {
+    if (cells(i/2, i%2) != INVALID) return;
+    std::queue<size_t> Q;
+    Q.push(i);
+    cells(i/2, i%2) = count;
+    while (!Q.empty()) {
+      const size_t index = Q.front();
       Q.pop();
-      const size_t patch_id = entry.first;
-      const short side = entry.second;
-      // face-edges adjacent to patch
-      const auto& adj_edges = patch_edge_adj[patch_id];
-      // Loop over **ALL EDGES IN THE ENTIRE PATCH** not even just the boundary
-      // edges, all edges... O(n)
-      for (const auto& ei : adj_edges)
-      {
-        // unique edge
-        const size_t uei = EMAP[ei];
-        // ordering of face-edges stored at edge
-        const auto& order = orders[uei];
-        // consistent orientation flags at face-edges stored at edge
-        const auto& orientation = orientations[uei];
-        const size_t edge_valance = order.size();
-        // Search for ei's (i.e. patch_id's) place in the ordering: O(#patches)
-        size_t curr_i = 0;
-        for (curr_i=0; curr_i < edge_valance; curr_i++)
-        {
-          if ((size_t)order[curr_i] == ei) break;
-        }
-        assert(curr_i < edge_valance && "Failed to find edge.");
-        // is the starting face consistent?
-        const bool cons = orientation[curr_i];
-        // Look clockwise or counter-clockwise for the next face, depending on
-        // whether looking to left or right side and whether consistently
-        // oriented or not.
-        size_t next;
-        if (side == 0)
-        {
-          next = (cons ? (curr_i + 1) :
-          (curr_i + edge_valance - 1)) % edge_valance;
-        } else {
-          next = (cons ? (curr_i+edge_valance-1) : (curr_i+1))%edge_valance;
-        }
-        // Determine face-edge index of next face
-        const size_t next_ei = order[next];
-        // Determine whether next is consistently oriented
-        const bool next_cons = orientation[next];
-        // Determine patch of next
-        const size_t next_patch_id = P[next_ei % num_faces];
-        // Determine which side of patch cell_idx is on, based on whether the
-        // consistency of next matches the consistency of this patch and which
-        // side of this patch we're on.
-        const short next_patch_side = (next_cons != cons) ?  side:abs(side-1);
-        // If that side of next's patch is not labeled, then label and add to
-        // queue
-        if (cells(next_patch_id, next_patch_side) == INVALID) 
-        {
-          Q.emplace(next_patch_id, next_patch_side);
-          cells(next_patch_id, next_patch_side) = cell_idx;
-        }else 
-        {
-          assert(
-            (size_t)cells(next_patch_id, next_patch_side) == cell_idx && 
-            "Encountered cell assignment inconsistency");
+      for (const auto j : equivalent_cells[index]) {
+        if (cells(j/2, j%2) == INVALID) {
+          cells(j/2, j%2) = count;
+          Q.push(j);
         }
       }
     }
+    count++;
   };
-
-  size_t count=0;
-  // Loop over all patches
-  for (size_t i=0; i<num_patches; i++)
-  {
-    // If left side of patch is still unlabeled, create a new cell and find all
-    // patches also on its boundary
-    if (cells(i, 0) == INVALID)
-    {
-      peel_cell_bd(i, 0, count,cells);
-      count++;
-    }
-    // Likewise for right side
-    if (cells(i, 1) == INVALID)
-    {
-      peel_cell_bd(i, 1, count,cells);
-      count++;
-    }
+  for (size_t i=0; i<num_patches; i++) {
+    extract_equivalent_cells(i*2);
+    extract_equivalent_cells(i*2+1);
   }
+
+  assert((cells.array() != INVALID).all());
   return count;
 }
 

include/igl/copyleft/cgal/order_facets_around_edge.cpp

@@ -274,124 +274,126 @@ void igl::copyleft::cgal::order_facets_around_edge(
 }
 
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedI>
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI>
 IGL_INLINE
 void igl::copyleft::cgal::order_facets_around_edge(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        size_t s,
-        size_t d, 
-        const std::vector<int>& adj_faces,
-        const Eigen::PlainObjectBase<DerivedV>& pivot_point,
-        Eigen::PlainObjectBase<DerivedI>& order)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  size_t s,
+  size_t d, 
+  const std::vector<int>& adj_faces,
+  const Eigen::PlainObjectBase<DerivedV>& pivot_point,
+  Eigen::PlainObjectBase<DerivedI>& order)
 {
+  assert(V.cols() == 3);
+  assert(F.cols() == 3);
+  assert(pivot_point.cols() == 3);
+  auto signed_index_to_index = [&](int signed_idx)
+  {
+      return abs(signed_idx) -1;
+  };
+  auto get_opposite_vertex_index = [&](size_t fid) -> typename DerivedF::Scalar
+  {
+      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;
+  };
 
-    assert(V.cols() == 3);
-    assert(F.cols() == 3);
-    assert(pivot_point.cols() == 3);
-    auto signed_index_to_index = [&](int signed_idx)
-    {
-        return abs(signed_idx) -1;
-    };
-    auto get_opposite_vertex_index = [&](size_t fid) -> typename DerivedF::Scalar
-    {
-        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 std::runtime_error(
-                    "Pivot point is collinear with the outer edge!");
-        }
+  {
+    // 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 std::runtime_error(
+                "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
+  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(adj_faces[i]) <
-            signed_index_to_index(adj_faces[j]);
-    };
-    std::vector<size_t> adj_order(N);
-    for (size_t i=0; i<N; i++) adj_order[i] = i;
-    std::sort(adj_order.begin(), adj_order.end(), comp);
+  // 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(adj_faces[i]) <
+      signed_index_to_index(adj_faces[j]);
+  };
+  std::vector<size_t> adj_order(N);
+  for (size_t i=0; i<N; i++) adj_order[i] = i;
+  std::sort(adj_order.begin(), adj_order.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[adj_order[i]]);
-        vertices.row(i) = V.row(get_opposite_vertex_index(fid));
-    }
-    vertices.row(N  ) = pivot_point;
-    vertices.row(N+1) = V.row(s);
-    vertices.row(N+2) = V.row(d);
+  DerivedV vertices(num_faces + 2, 3);
+  for (size_t i=0; i<N; i++) 
+  {
+    const size_t fid = signed_index_to_index(adj_faces[adj_order[i]]);
+    vertices.row(i) = V.row(get_opposite_vertex_index(fid));
+  }
+  vertices.row(N  ) = pivot_point;
+  vertices.row(N+1) = V.row(s);
+  vertices.row(N+2) = V.row(d);
 
-    DerivedF faces(num_faces, 3);
-    for (size_t i=0; i<N; i++)
+  DerivedF faces(num_faces, 3);
+  for (size_t i=0; i<N; i++)
+  {
+    if (adj_faces[adj_order[i]] < 0) 
     {
-        if (adj_faces[adj_order[i]] < 0) {
-            faces(i,0) = N+1; // s
-            faces(i,1) = N+2; // d
-            faces(i,2) = i  ;
-        } else {
-            faces(i,0) = N+2; // d
-            faces(i,1) = N+1; // s
-            faces(i,2) = i  ;
-        }
+      faces(i,0) = N+1; // s
+      faces(i,1) = N+2; // d
+      faces(i,2) = i  ;
+    } else 
+    {
+      faces(i,0) = N+2; // d
+      faces(i,1) = N+1; // s
+      faces(i,2) = i  ;
     }
-    // Last face is the pivot face.
-    faces(N, 0) = N+1;
-    faces(N, 1) = N+2;
-    faces(N, 2) = N;
+  }
+  // Last face is the pivot face.
+  faces(N, 0) = N+1;
+  faces(N, 1) = N+2;
+  faces(N, 2) = N;
 
-    std::vector<int> adj_faces_with_pivot(num_faces);
-    for (size_t i=0; i<num_faces; i++)
+  std::vector<int> adj_faces_with_pivot(num_faces);
+  for (size_t i=0; i<num_faces; i++)
+  {
+    if ((size_t)faces(i,0) == N+1 && (size_t)faces(i,1) == N+2)
     {
-        if ((size_t)faces(i,0) == N+1 && (size_t)faces(i,1) == N+2)
-        {
-            adj_faces_with_pivot[i] = int(i+1) * -1;
-        } else
-        {
-            adj_faces_with_pivot[i] = int(i+1);
-        }
+        adj_faces_with_pivot[i] = int(i+1) * -1;
+    } else
+    {
+        adj_faces_with_pivot[i] = int(i+1);
     }
+  }
 
-    DerivedI order_with_pivot;
-    order_facets_around_edge(
-            vertices, faces, N+1, N+2,
-            adj_faces_with_pivot, order_with_pivot);
+  DerivedI order_with_pivot;
+  order_facets_around_edge(
+    vertices, faces, N+1, N+2, adj_faces_with_pivot, order_with_pivot);
 
-    assert((size_t)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++)
+  assert((size_t)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 ((size_t)order_with_pivot[i] == N)
     {
-        if ((size_t)order_with_pivot[i] == N)
-        {
-            pivot_index = i;
-            break;
-        }
+      pivot_index = i;
+      break;
     }
-    assert(pivot_index < num_faces);
+  }
+  assert(pivot_index < num_faces);
 
-    for (size_t i=0; i<N; i++)
-    {
-        order[i] = adj_order[order_with_pivot[(pivot_index+i+1)%num_faces]];
-    }
+  for (size_t i=0; i<N; i++)
+  {
+    order[i] = adj_order[order_with_pivot[(pivot_index+i+1)%num_faces]];
+  }
 }
 
 

include/igl/copyleft/cgal/order_facets_around_edge.h

@@ -12,10 +12,12 @@
 #include <Eigen/Core>
 #include <vector>
 
-namespace igl {
+namespace igl 
+{
   namespace copyleft
   {
-    namespace cgal {
+    namespace cgal 
+    {
       // Given a directed edge, sort its adjacent faces.  Assuming the
       // directed edge is (s, d).  Sort the adjacent faces clockwise around the
       // axis (d - s), i.e. left-hand rule.  An adjacent face is consistently
@@ -26,15 +28,14 @@ namespace igl {
       // computed as (consistent? 1:-1) * (face_index + 1).
       //
       // Inputs:
-      //   V          #V by 3 list of vertices.
-      //   F          #F by 3 list of faces
-      //   s          Index of source vertex.
-      //   d          Index of desination vertex.
+      //   V  #V by 3 list of vertices.
+      //   F  #F by 3 list of faces
+      //   s  Index of source vertex.
+      //   d  Index of desination vertex.
       //   adj_faces  List of adjacent face signed indices.
-      //
       // Output:
-      //   order      List of face indices that orders adjacent faces around
-      //              edge (s, d) clockwise.
+      //   order  List of face indices that orders adjacent faces around edge
+      //     (s, d) clockwise.
       template<
         typename DerivedV,
         typename DerivedF,
@@ -43,7 +44,8 @@ namespace igl {
       void 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,
           bool debug=false);
@@ -58,12 +60,13 @@ namespace igl {
         typename DerivedI>
       IGL_INLINE
       void order_facets_around_edge(
-          const Eigen::PlainObjectBase<DerivedV>& V,
-          const Eigen::PlainObjectBase<DerivedF>& F,
-          size_t s, size_t d, 
-          const std::vector<int>& adj_faces,
-          const Eigen::PlainObjectBase<DerivedV>& pivot_point,
-          Eigen::PlainObjectBase<DerivedI>& order);
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        size_t s, 
+        size_t d, 
+        const std::vector<int>& adj_faces,
+        const Eigen::PlainObjectBase<DerivedV>& pivot_point,
+        Eigen::PlainObjectBase<DerivedI>& order);
     }
   }
 }

include/igl/copyleft/cgal/outer_hull.cpp

@@ -5,8 +5,107 @@
 // 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 "points_inside_component.h"
 #include "outer_hull.h"
+#include "extract_cells.h"
+#include "remesh_self_intersections.h"
+#include "../../remove_unreferenced.h"
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedHV,
+  typename DerivedHF,
+  typename DerivedJ,
+  typename Derivedflip>
+IGL_INLINE void igl::copyleft::cgal::outer_hull(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedHV> & HV,
+  Eigen::PlainObjectBase<DerivedHF> & HF,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<Derivedflip> & flip)
+{
+  // Exact types
+  typedef CGAL::Epeck Kernel;
+  typedef Kernel::FT ExactScalar;
+  typedef 
+    Eigen::Matrix<
+    ExactScalar,
+    Eigen::Dynamic,
+    Eigen::Dynamic,
+    DerivedHV::IsRowMajor> 
+      MatrixXES;
+  // Remesh self-intersections
+  MatrixXES Vr;
+  DerivedHF Fr;
+  DerivedJ Jr;
+  {
+    RemeshSelfIntersectionsParam params;
+    Eigen::VectorXi I;
+    Eigen::MatrixXi IF;
+    remesh_self_intersections(V, F, params, Vr, Fr, IF, Jr, I);
+    // Merge coinciding vertices into non-manifold vertices.
+    std::for_each(Fr.data(), Fr.data()+Fr.size(),
+      [&I](typename DerivedHF::Scalar& a) { a=I[a]; });
+      // Remove unreferenced vertices.
+    Eigen::VectorXi UIM;
+    remove_unreferenced(MatrixXES(Vr),DerivedHF(Fr), Vr, Fr, UIM);
+  }
+  // Extract cells for each face
+  Eigen::MatrixXi C;
+  extract_cells(Vr,Fr,C);
+  // Extract faces on ambient cell
+  int num_outer = 0;
+  for(int i = 0;i<C.rows();i++)
+  {
+    num_outer += ( C(i,0) == 0 || C(i,1) == 0 ) ? 1 : 0;
+  }
+  HF.resize(num_outer,3);
+  J.resize(num_outer,1);
+  flip.resize(num_outer,1);
+  {
+    int h = 0;
+    for(int i = 0;i<C.rows();i++)
+    {
+      if(C(i,0)==0)
+      {
+        HF.row(h) = Fr.row(i);
+        J(h) = Jr(i);
+        flip(h) = false;
+        h++;
+      }else if(C(i,1) == 0)
+      {
+        HF.row(h) = Fr.row(i).reverse();
+        J(h) = Jr(i);
+        flip(h) = true;
+        h++;
+      }
+    }
+    assert(h == num_outer);
+  }
+  // Remove unreferenced vertices and re-index faces
+  {
+    // Cast to output type
+    DerivedHV Vr_cast(Vr.rows(),Vr.cols());
+    for(int i = 0;i<Vr.rows();i++)
+    {
+      for(int j = 0;j<Vr.cols();j++)
+      {
+        assign_scalar(Vr(i,j), Vr_cast(i,j));
+      }
+    }
+    Eigen::VectorXi I;
+    remove_unreferenced(Vr_cast,DerivedHF(HF),HV,HF,I);
+  }
+}
+
+#include "points_inside_component.h"
 #include "order_facets_around_edges.h"
 #include "outer_facet.h"
 #include "../../sortrows.h"
@@ -16,10 +115,7 @@
 #include "../../unique_edge_map.h"
 #include "../../barycenter.h"
 #include "../../per_face_normals.h"
-#include "../../writePLY.h"
 #include "../../sort_angles.h"
-#include "../../writePLY.h"
-
 #include <Eigen/Geometry>
 #include <vector>
 #include <map>
@@ -35,7 +131,7 @@ template <
   typename DerivedG,
   typename DerivedJ,
   typename Derivedflip>
-IGL_INLINE void igl::copyleft::cgal::outer_hull(
+IGL_INLINE void igl::copyleft::cgal::outer_hull_legacy(
   const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::PlainObjectBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedG> & G,
@@ -435,6 +531,4 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::copyleft::cgal::outer_hull<Eigen::Matrix<double, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-template void igl::copyleft::cgal::outer_hull<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::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/outer_hull.h

@@ -15,6 +15,32 @@ namespace igl
   {
     namespace cgal
     {
+      // Compute the "outer hull" of a piecewise constant winding number induce
+      // triangle mesh (V,F).
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions
+      //   F  #F by 3 list of triangle indices into V
+      // Outputs:
+      //   HV  #HV by 3 list of output vertex positions
+      //   HF  #HF by 3 list of output triangle indices into HV
+      //   J  #HF list of indices into F
+      //   flip  #HF list of whether facet was flipped when added to HF
+      //
+      template <
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedHV,
+        typename DerivedHF,
+        typename DerivedJ,
+        typename Derivedflip>
+      IGL_INLINE void outer_hull(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<DerivedHV> & HV,
+        Eigen::PlainObjectBase<DerivedHF> & HF,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<Derivedflip> & flip);
       // Compute the "outer hull" of a potentially non-manifold mesh (V,F) whose
       // intersections have been "resolved" (e.g. using `cork` or
       // `igl::copyleft::cgal::selfintersect`). The outer hull is defined to be all facets
@@ -24,6 +50,10 @@ namespace igl
       // "flaps".  This implementation largely follows Section 3.6 of "Direct
       // repair of self-intersecting meshes" [Attene 2014].
       //
+      // Note: This doesn't require the input mesh to be piecewise constant
+      // winding number, but won't handle multiple non-nested connected
+      // components.
+      //
       // Inputs:
       //   V  #V by 3 list of vertex positions
       //   F  #F by 3 list of triangle indices into V
@@ -38,7 +68,7 @@ namespace igl
         typename DerivedG,
         typename DerivedJ,
         typename Derivedflip>
-      IGL_INLINE void outer_hull(
+      IGL_INLINE void outer_hull_legacy(
         const Eigen::PlainObjectBase<DerivedV> & V,
         const Eigen::PlainObjectBase<DerivedF> & F,
         Eigen::PlainObjectBase<DerivedG> & G,

include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -74,7 +74,7 @@ IGL_INLINE size_t igl::copyleft::cgal::peel_outer_hull_layers(
       writePLY(ss.str(), vertices, Fr);
   }
 #endif
-    outer_hull(V,Fr,Fo,Jo,flipr);
+    outer_hull_legacy(V,Fr,Fo,Jo,flipr);
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
   writePLY(STR("outer-hull-output-"<<iter<<".ply"),V,Fo);
   cout<<"reindex, flip..."<<endl;

include/igl/copyleft/cgal/piecewise_constant_winding_number.cpp

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

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -34,11 +34,12 @@ template<
   typename DerivedF,
   typename DerivedL,
   typename DerivedW>
-IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
+IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     const Eigen::PlainObjectBase<DerivedL>& labels,
-    Eigen::PlainObjectBase<DerivedW>& W) {
+    Eigen::PlainObjectBase<DerivedW>& W) 
+{
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
   const auto & tictoc = []() -> double
   {
@@ -60,9 +61,11 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
   Eigen::VectorXi EMAP;
   std::vector<std::vector<size_t> > uE2E;
   igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  bool valid = true;
   if (!piecewise_constant_winding_number(F, uE, uE2E)) 
   {
     std::cerr << "Input mesh is not orientable!" << std::endl;
+    valid = false;
   }
 
   Eigen::VectorXi P;
@@ -138,26 +141,29 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
           Q.pop();
           int curr_label = cell_labels[curr_idx];
           for (const auto& neighbor : cell_adjacency[curr_idx]) {
-            if (cell_labels[std::get<0>(neighbor)] == 0) {
+            if (cell_labels[std::get<0>(neighbor)] == 0) 
+            {
               cell_labels[std::get<0>(neighbor)] = curr_label * -1;
               Q.push(std::get<0>(neighbor));
               parents[std::get<0>(neighbor)] = curr_idx;
-            } else {
-              if (cell_labels[std::get<0>(neighbor)] !=
-                  curr_label * -1) {
+            } else 
+            {
+              if (cell_labels[std::get<0>(neighbor)] != curr_label * -1) 
+              {
                 std::cerr << "Odd cell cycle detected!" << std::endl;
                 auto path = trace_parents(curr_idx);
                 path.reverse();
                 auto path2 = trace_parents(std::get<0>(neighbor));
-                path.insert(path.end(),
-                    path2.begin(), path2.end());
-                for (auto cell_id : path) {
+                path.insert(path.end(), path2.begin(), path2.end());
+                for (auto cell_id : path) 
+                {
                   std::cout << cell_id << " ";
                   std::stringstream filename;
                   filename << "cell_" << cell_id << ".ply";
                   save_cell(filename.str(), cell_id);
                 }
                 std::cout << std::endl;
+                valid = false;
               }
               // Do not fail when odd cycle is detected because the resulting
               // integer winding number field, although inconsistent, may still
@@ -249,7 +255,8 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
 #endif
 
   W.resize(num_faces, num_labels*2);
-  for (size_t i=0; i<num_faces; i++) {
+  for (size_t i=0; i<num_faces; i++) 
+  {
     const size_t patch = P[i];
     const size_t positive_cell = per_patch_cells(patch, 0);
     const size_t negative_cell = per_patch_cells(patch, 1);
@@ -261,9 +268,10 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
   log_time("store_result");
 #endif
+  return valid;
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::copyleft::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<int, -1, -1, 0, -1, -1> >&);
+template bool igl::copyleft::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/copyleft/cgal/propagate_winding_numbers.h

@@ -19,11 +19,16 @@
 // clockwise facing facets should equal the number of counterclockwise facing
 // facets.
 
-namespace igl {
+namespace igl
+{
   namespace copyleft
   {
-    namespace cgal {
-
+    namespace cgal
+    {
+      // TODO: This shouldn't need to be in igl::copyleft::cgal, it should
+      // instead take as input an index of the ambient cell and the winding
+      // number vector there.
+      //
       // Compute winding number on each side of the face.  The input mesh
       // could contain multiple connected components.  The input mesh must
       // represent the boundary of a valid 3D volume, which means it is
@@ -33,23 +38,24 @@ namespace igl {
       //   V  #V by 3 list of vertex positions.
       //   F  #F by 3 list of triangle indices into V.
       //   labels  #F list of facet labels ranging from 0 to k-1.
-      //
       // Output:
       //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
       //      number on the positive side of facet ``i`` with respect to the
       //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
       //      number on the negative side of facet ``i`` with respect to the
       //      facets labeled ``j``.
+      // Returns true iff the input induces a piecewise-constant winding number
+      //   field.
       template<
         typename DerivedV,
         typename DerivedF,
         typename DerivedL,
         typename DerivedW>
-      IGL_INLINE void propagate_winding_numbers(
-          const Eigen::PlainObjectBase<DerivedV>& V,
-          const Eigen::PlainObjectBase<DerivedF>& F,
-          const Eigen::PlainObjectBase<DerivedL>& labels,
-          Eigen::PlainObjectBase<DerivedW>& W);
+      IGL_INLINE bool propagate_winding_numbers(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedL>& labels,
+        Eigen::PlainObjectBase<DerivedW>& W);
     }
   }
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -19,6 +19,33 @@
 
 //#define REMESH_INTERSECTIONS_TIMING
 
+template <
+typename DerivedV,
+typename DerivedF,
+typename Kernel,
+typename DerivedVV,
+typename DerivedFF,
+typename DerivedJ,
+typename DerivedIM>
+IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const std::vector<CGAL::Triangle_3<Kernel> > & T,
+        const std::map<
+        typename DerivedF::Index,
+        std::vector<
+        std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
+        const std::map<
+        std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+        std::vector<typename DerivedF::Index> > & edge2faces,
+        Eigen::PlainObjectBase<DerivedVV> & VV,
+        Eigen::PlainObjectBase<DerivedFF> & FF,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM) {
+    igl::copyleft::cgal::remesh_intersections(
+            V, F, T, offending, edge2faces, false, VV, FF, J, IM);
+}
+
 template <
 typename DerivedV,
 typename DerivedF,
@@ -38,6 +65,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         const std::map<
         std::pair<typename DerivedF::Index,typename DerivedF::Index>,
         std::vector<typename DerivedF::Index> > & /*edge2faces*/,
+        bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,
@@ -187,51 +215,61 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
      * Given p on triangle indexed by ori_f, determine the index of p.
      */
     auto determine_point_index = [&](
-            const Point_3& p, const size_t ori_f) -> Index {
+        const Point_3& p, const size_t ori_f) -> Index {
+      if (stitch_all) {
+        // No need to check if p shared by multiple triangles because all shared
+        // vertices would be merged later on.
+        const size_t index = num_base_vertices + new_vertices.size();
+        new_vertices.push_back(p);
+        return index;
+      } else {
+        // Stitching triangles according to input connectivity.
+        // This step is potentially costly.
         const auto& triangle = T[ori_f];
         const auto& f = F.row(ori_f).eval();
 
         // Check if p is one of the triangle corners.
         for (size_t i=0; i<3; i++) {
-            if (p == triangle[i]) return f[i];
+          if (p == triangle[i]) return f[i];
         }
 
         // Check if p is on one of the edges.
         for (size_t i=0; i<3; i++) {
-            const Point_3 curr_corner = triangle[i];
-            const Point_3 next_corner = triangle[(i+1)%3];
-            const Segment_3 edge(curr_corner, next_corner);
-            if (edge.has_on(p)) {
-                const Index curr = f[i];
-                const Index next = f[(i+1)%3];
-                Edge key;
-                key.first = curr<next?curr:next;
-                key.second = curr<next?next:curr;
-                auto itr = edge_vertices.find(key);
-                if (itr == edge_vertices.end()) {
-                    const Index index =
-                        num_base_vertices + new_vertices.size();
-                    edge_vertices.insert({key, {index}});
-                    new_vertices.push_back(p);
-                    return index;
-                } else {
-                    for (const auto vid : itr->second) {
-                        if (p == new_vertices[vid - num_base_vertices]) {
-                            return vid;
-                        }
-                    }
-                    const size_t index = num_base_vertices + new_vertices.size();
-                    new_vertices.push_back(p);
-                    itr->second.push_back(index);
-                    return index;
+          const Point_3 curr_corner = triangle[i];
+          const Point_3 next_corner = triangle[(i+1)%3];
+          const Segment_3 edge(curr_corner, next_corner);
+          if (edge.has_on(p)) {
+            const Index curr = f[i];
+            const Index next = f[(i+1)%3];
+            Edge key;
+            key.first = curr<next?curr:next;
+            key.second = curr<next?next:curr;
+            auto itr = edge_vertices.find(key);
+            if (itr == edge_vertices.end()) {
+              const Index index =
+                num_base_vertices + new_vertices.size();
+              edge_vertices.insert({key, {index}});
+              new_vertices.push_back(p);
+              return index;
+            } else {
+              for (const auto vid : itr->second) {
+                if (p == new_vertices[vid - num_base_vertices]) {
+                  return vid;
                 }
+              }
+              const size_t index = num_base_vertices + new_vertices.size();
+              new_vertices.push_back(p);
+              itr->second.push_back(index);
+              return index;
             }
+          }
         }
 
         // p must be in the middle of the triangle.
         const size_t index = num_base_vertices + new_vertices.size();
         new_vertices.push_back(p);
         return index;
+      }
     };
 
     /**
@@ -241,32 +279,45 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             const std::vector<Point_3>& vertices,
             const std::vector<std::vector<Index> >& faces,
             const std::vector<Index>& involved_faces) -> void {
-        for (const auto& f : faces) {
-            const Point_3& v0 = vertices[f[0]];
-            const Point_3& v1 = vertices[f[1]];
-            const Point_3& v2 = vertices[f[2]];
-            Point_3 center(
-                    (v0[0] + v1[0] + v2[0]) / 3.0,
-                    (v0[1] + v1[1] + v2[1]) / 3.0,
-                    (v0[2] + v1[2] + v2[2]) / 3.0);
-            for (const auto& ori_f : involved_faces) {
-                const auto& triangle = T[ori_f];
-                const Plane_3 P = triangle.supporting_plane();
-                if (triangle.has_on(center)) {
-                    std::vector<Index> corners(3);
-                    corners[0] = determine_point_index(v0, ori_f);
-                    corners[1] = determine_point_index(v1, ori_f);
-                    corners[2] = determine_point_index(v2, ori_f);
-                    if (CGAL::orientation(
-                                P.to_2d(v0), P.to_2d(v1), P.to_2d(v2))
-                            == CGAL::RIGHT_TURN) {
-                        std::swap(corners[0], corners[1]);
-                    }
-                    resolved_faces.emplace_back(corners);
-                    source_faces.push_back(ori_f);
-                }
+      assert(involved_faces.size() > 0);
+      for (const auto& f : faces) {
+        const Point_3& v0 = vertices[f[0]];
+        const Point_3& v1 = vertices[f[1]];
+        const Point_3& v2 = vertices[f[2]];
+        Point_3 center(
+            (v0[0] + v1[0] + v2[0]) / 3.0,
+            (v0[1] + v1[1] + v2[1]) / 3.0,
+            (v0[2] + v1[2] + v2[2]) / 3.0);
+        if (involved_faces.size() == 1) {
+          // If only there is only one involved face, all sub-triangles must
+          // belong to it and have the correct orientation.
+          const auto& ori_f = involved_faces[0];
+          std::vector<Index> corners(3);
+          corners[0] = determine_point_index(v0, ori_f);
+          corners[1] = determine_point_index(v1, ori_f);
+          corners[2] = determine_point_index(v2, ori_f);
+          resolved_faces.emplace_back(corners);
+          source_faces.push_back(ori_f);
+        } else {
+          for (const auto& ori_f : involved_faces) {
+            const auto& triangle = T[ori_f];
+            const Plane_3 P = triangle.supporting_plane();
+            if (triangle.has_on(center)) {
+              std::vector<Index> corners(3);
+              corners[0] = determine_point_index(v0, ori_f);
+              corners[1] = determine_point_index(v1, ori_f);
+              corners[2] = determine_point_index(v2, ori_f);
+              if (CGAL::orientation(
+                    P.to_2d(v0), P.to_2d(v1), P.to_2d(v2))
+                  == CGAL::RIGHT_TURN) {
+                std::swap(corners[0], corners[1]);
+              }
+              resolved_faces.emplace_back(corners);
+              source_faces.push_back(ori_f);
             }
+          }
         }
+      }
     };
 
     // Process un-touched faces.
@@ -335,12 +386,11 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     log_time("cdt");
 #endif
 
-    // Post process
     for (size_t i=0; i<num_cdts; i++) {
-        const auto& vertices = cdt_vertices[i];
-        const auto& faces = cdt_faces[i];
-        const auto& involved_faces = cdt_inputs[i].second;
-        post_triangulation_process(vertices, faces, involved_faces);
+      const auto& vertices = cdt_vertices[i];
+      const auto& faces = cdt_faces[i];
+      const auto& involved_faces = cdt_inputs[i].second;
+      post_triangulation_process(vertices, faces, involved_faces);
     }
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("stitching");
@@ -371,28 +421,27 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     J.resize(num_out_faces);
     std::copy(source_faces.begin(), source_faces.end(), J.data());
 
-    struct PointHash {
-        size_t operator()(const Point_3& p) const {
-            static const std::hash<double> hasher{};
-            double x,y,z;
-            assign_scalar(p.x(), x);
-            assign_scalar(p.y(), y);
-            assign_scalar(p.z(), z);
-            return hasher(x) ^ hasher(y) ^ hasher(z);
-        }
-    };
-
-    // TODO: use igl::unique()
     // Extract unique vertex indices.
-
     const size_t VV_size = VV.rows();
     IM.resize(VV_size,1);
 
     DerivedVV unique_vv;
     Eigen::VectorXi unique_to_vv, vv_to_unique;
     igl::unique_rows(VV, unique_vv, unique_to_vv, vv_to_unique);
-    for (Index v=0; v<VV_size; v++) {
-      IM(v) = unique_to_vv[vv_to_unique[v]];
+    if (!stitch_all) {
+      // Vertices with the same coordinates would be represented by one vertex.
+      // The IM value of an vertex is the index of the representative vertex.
+      for (Index v=0; v<VV_size; v++) {
+        IM(v) = unique_to_vv[vv_to_unique[v]];
+      }
+    } else {
+      // Screw IM and representative vertices.  Merge all vertices having the
+      // same coordinates into a single vertex and set IM to identity map.
+      VV = unique_vv;
+      std::transform(FF.data(), FF.data() + FF.rows()*FF.cols(),
+          FF.data(), [&vv_to_unique](const typename DerivedFF::Scalar& a)
+          { return vv_to_unique[a]; });
+      IM.setLinSpaced(unique_vv.rows(), 0, unique_vv.rows()-1);
     }
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("store_results");

include/igl/copyleft/cgal/remesh_intersections.h

@@ -62,6 +62,35 @@ namespace igl
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,
         Eigen::PlainObjectBase<DerivedIM> & IM);
+
+      // Same as above except ``stitch_all`` flag:
+      //
+      // Input:
+      //   stitch_all: if true, merge all vertices with thte same coordiante.
+      template <
+        typename DerivedV,
+        typename DerivedF,
+        typename Kernel,
+        typename DerivedVV,
+        typename DerivedFF,
+        typename DerivedJ,
+        typename DerivedIM>
+      IGL_INLINE void remesh_intersections(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const std::vector<CGAL::Triangle_3<Kernel> > & T,
+        const std::map<
+          typename DerivedF::Index,
+            std::vector<
+            std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
+        const std::map<
+          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+          std::vector<typename DerivedF::Index> > & edge2faces,
+        bool stitch_all,
+        Eigen::PlainObjectBase<DerivedVV> & VV,
+        Eigen::PlainObjectBase<DerivedFF> & FF,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM);
     }
   }
 }

include/igl/copyleft/cgal/submesh_aabb_tree.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// 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 "submesh_aabb_tree.h"
+#include <stdexcept>
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename Kernel>
+IGL_INLINE void igl::copyleft::cgal::submesh_aabb_tree(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedI>& I,
+  CGAL::AABB_tree<
+    CGAL::AABB_traits<
+      Kernel, 
+      CGAL::AABB_triangle_primitive<
+        Kernel, typename std::vector<
+          typename Kernel::Triangle_3 >::iterator > > > & tree,
+  std::vector<typename Kernel::Triangle_3 > & triangles,
+  std::vector<bool> & in_I)
+{
+  in_I.resize(F.rows(), false);
+  const size_t num_faces = I.rows();
+  for (size_t i=0; i<num_faces; i++) 
+  {
+    const Eigen::Vector3i f = F.row(I(i, 0));
+    in_I[I(i,0)] = true;
+    triangles.emplace_back(
+      typename Kernel::Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+      typename Kernel::Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+      typename Kernel::Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+#ifndef NDEBUG
+    if (triangles.back().is_degenerate()) 
+    {
+      throw std::runtime_error(
+        "Input facet components contains degenerated triangles");
+    }
+#endif
+  }
+  tree.insert(triangles.begin(), triangles.end());
+  tree.accelerate_distance_queries();
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::copyleft::cgal::submesh_aabb_tree<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epeck, CGAL::AABB_triangle_primitive<CGAL::Epeck, std::vector<CGAL::Epeck::Triangle_3, std::allocator<CGAL::Epeck::Triangle_3> >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Epeck::Triangle_3, std::allocator<CGAL::Epeck::Triangle_3> >&, std::vector<bool, std::allocator<bool> >&);
+#endif

include/igl/copyleft/cgal/submesh_aabb_tree.h

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// 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_COPYLET_CGAL_SUBMESH_AABB_TREE_H
+#define IGL_COPYLET_CGAL_SUBMESH_AABB_TREE_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Build an AABB tree for a submesh indicated by a face selection list I
+      // of a full mesh (V,F)
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertices.
+      //   F  #F by 3 array of faces.
+      //   I  #I list of triangle indices to consider.
+      // Outputs:
+      //   tree  aabb containing triangles of (V,F(I,:))
+      //   triangles  #I list of cgal triangles
+      //   in_I  #F list of whether in submesh
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename Kernel>
+      IGL_INLINE void submesh_aabb_tree(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        CGAL::AABB_tree<
+          CGAL::AABB_traits<
+            Kernel, 
+            CGAL::AABB_triangle_primitive<
+              Kernel, typename std::vector<
+                typename Kernel::Triangle_3 >::iterator > > > & tree,
+        std::vector<typename Kernel::Triangle_3 > & triangles,
+        std::vector<bool> & in_I);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "submesh_aabb_tree.cpp"
+#endif
+
+#endif

include/igl/edge_lengths.cpp

@@ -99,4 +99,5 @@ template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::M
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/embree/unproject_onto_mesh.cpp

@@ -7,8 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "unproject_onto_mesh.h"
 #include "EmbreeIntersector.h"
-#include <igl/unproject.h>
-#include <igl/embree/unproject_in_mesh.h>
+#include "../unproject_onto_mesh.h"
 #include <vector>
 
 IGL_INLINE bool igl::embree::unproject_onto_mesh(
@@ -23,20 +22,14 @@ IGL_INLINE bool igl::embree::unproject_onto_mesh(
 {
   using namespace std;
   using namespace Eigen;
-  MatrixXd obj;
-  vector<igl::Hit> hits;
-
-  // This is lazy, it will find more than just the first hit
-  unproject_in_mesh(pos,model,proj,viewport,ei,obj,hits);
-
-  if (hits.size()> 0)
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    igl::Hit & hit)->bool
   {
-    bc << 1.0-hits[0].u-hits[0].v, hits[0].u, hits[0].v;
-    fid = hits[0].id;
-    return true;
-  }
-
-  return false;
+    return ei.intersectRay(s,dir,hit);
+  };
+  return igl::unproject_onto_mesh(pos,model,proj,viewport,shoot_ray,fid,bc);
 }
 
 IGL_INLINE bool igl::embree::unproject_onto_mesh(
@@ -50,14 +43,14 @@ IGL_INLINE bool igl::embree::unproject_onto_mesh(
   int& vid)
 {
   Eigen::Vector3f bc;
-  bool hit = unproject_onto_mesh(pos,F,model,proj,viewport,ei,fid,bc);
-  int i;
-  if (hit)
+  if(!igl::embree::unproject_onto_mesh(pos,F,model,proj,viewport,ei,fid,bc))
   {
-    bc.maxCoeff(&i);
-    vid = F(fid,i);
+    return false;
   }
-  return hit;
+  int i;
+  bc.maxCoeff(&i);
+  vid = F(fid,i);
+  return true;
 }
 
 

include/igl/extract_manifold_patches.h

@@ -6,13 +6,9 @@
 #include <vector>
 
 namespace igl {
-    // Extract a set of maximal manifold patches from a given mesh.
-    // A maximal manifold patch is a subset of the input faces that is
-    // connected and manifold, and it cannot be expanded further by
-    // including more faces.
-    //
-    // Assumes the input mesh have all self-intersection resolved.  See
-    // ``igl::cgal::remesh_self_intersection`` for more details.
+    // Extract a set of maximal patches from a given mesh.
+    // A maximal patch is a subset of the input faces that are connected via
+    // manifold edges; a patch is as large as possible.
     //
     // Inputs:
     //   F  #F by 3 list representing triangles.
@@ -25,15 +21,15 @@ namespace igl {
     // Returns:
     //   number of manifold patches.
     template <
-        typename DerivedF,
-        typename DerivedEMAP,
-        typename uE2EType,
-        typename DerivedP>
+      typename DerivedF,
+      typename DerivedEMAP,
+      typename uE2EType,
+      typename DerivedP>
     IGL_INLINE size_t extract_manifold_patches(
-            const Eigen::PlainObjectBase<DerivedF>& F,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-            const std::vector<std::vector<uE2EType> >& uE2E,
-            Eigen::PlainObjectBase<DerivedP>& P);
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+      const std::vector<std::vector<uE2EType> >& uE2E,
+      Eigen::PlainObjectBase<DerivedP>& P);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "extract_manifold_patches.cpp"

include/igl/ply.h.REMOVED.git-id

@@ -1 +1 @@
-a11a359439a307b36f2320e726048ce3dd344fda
+aa5b1239a610dbd02506563452ec69172ae4165a

include/igl/point_simplex_squared_distance.cpp

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

include/igl/point_simplex_squared_distance.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_POINT_SIMPLEX_SQUARED_DISTANCE_H
+#define IGL_POINT_SIMPLEX_SQUARED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Determine squared distance from a point to linear simplex
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::PlainObjectBase<Derivedp> & p,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_simplex_squared_distance.cpp"
+#endif
+#endif

include/igl/pseudonormal_test.cpp

@@ -6,7 +6,8 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "pseudonormal_test.h"
-#include "AABB.h"
+#include "barycentric_coordinates.h"
+#include "doublearea.h"
 #include "project_to_line_segment.h"
 #include <cassert>
 
@@ -49,7 +50,7 @@ IGL_INLINE void igl::pseudonormal_test(
   const double epsilon = 1e-12;
   if(area>MIN_DOUBLE_AREA)
   {
-    AABB<Eigen::MatrixXd,3>::barycentric_coordinates( c,A,B,C,b);
+    barycentric_coordinates( c,A,B,C,b);
     // Determine which normal to use
     const int type = (b.array()<=epsilon).cast<int>().sum();
     switch(type)

include/igl/ray_mesh_intersect.cpp

@@ -11,10 +11,10 @@ template <
   typename DerivedV, 
   typename DerivedF> 
 IGL_INLINE bool igl::ray_mesh_intersect(
-  const Eigen::PlainObjectBase<Derivedsource> & s,
-  const Eigen::PlainObjectBase<Deriveddir> & dir,
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedsource> & s,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   std::vector<igl::Hit> & hits)
 {
   using namespace Eigen;
@@ -53,10 +53,10 @@ template <
   typename DerivedV, 
   typename DerivedF> 
 IGL_INLINE bool igl::ray_mesh_intersect(
-  const Eigen::PlainObjectBase<Derivedsource> & source,
-  const Eigen::PlainObjectBase<Deriveddir> & dir,
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedsource> & source,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   igl::Hit & hit)
 {
   std::vector<igl::Hit> hits;
@@ -70,3 +70,7 @@ IGL_INLINE bool igl::ray_mesh_intersect(
     return false;
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::__1::vector<igl::Hit, std::__1::allocator<igl::Hit> >&);
+#endif

include/igl/ray_mesh_intersect.h

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

include/igl/readPLY.cpp

@@ -153,7 +153,6 @@ IGL_INLINE bool igl::readPLY(
     }
   }
   ply_close(in_ply);
-  fclose(fp);
   return true;
 }
 

include/igl/rotation_matrix_from_directions.cpp

@@ -11,20 +11,14 @@
 #include <iostream>
 
 template <typename Scalar>
-IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(const Eigen::Matrix<Scalar, 3, 1> v0,
-                                                                        const Eigen::Matrix<Scalar, 3, 1> v1,
-                                                                        bool normalized)
+IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(
+  const Eigen::Matrix<Scalar, 3, 1> v0,
+  const Eigen::Matrix<Scalar, 3, 1> v1)
 {
   Eigen::Matrix<Scalar, 3, 3> rotM;
   const double epsilon=1e-8;
-  // if (!normalized)
-//   {
-    // v0.normalize();
-    // v1.normalize();
-  // }
   Scalar dot=v0.normalized().dot(v1.normalized());
   ///control if there is no rotation
-//  if (dot>((double)1-epsilon))
   if ((v0-v1).norm()<epsilon)
   {
     rotM = Eigen::Matrix<Scalar, 3, 3>::Identity();
@@ -65,5 +59,5 @@ IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(cons
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template Eigen::Matrix<double, 3, 3, 0, 3, 3> igl::rotation_matrix_from_directions<double>(const Eigen::Matrix<double, 3, 1, 0, 3, 1>, const Eigen::Matrix<double, 3, 1, 0, 3, 1>, const bool);
+template Eigen::Matrix<double, 3, 3, 0, 3, 3> igl::rotation_matrix_from_directions<double>(const Eigen::Matrix<double, 3, 1, 0, 3, 1>, const Eigen::Matrix<double, 3, 1, 0, 3, 1>);
 #endif

include/igl/rotation_matrix_from_directions.h

@@ -1,36 +1,34 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>, Olga Diamanti <olga.diam@gmail.com>
+// Copyright (C) 2016 Alec Jacobson, Daniele Panozzo, Olga Diamanti 
 //
 // 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_ROTATION_MATRIX_FROM_DIRECTIONS
-#define IGL_ROTATION_MATRIX_FROM_DIRECTIONS
+#ifndef IGL_ROTATION_MATRIX_FROM_DIRECTIONS_H
+#define IGL_ROTATION_MATRIX_FROM_DIRECTIONS_H
 #include "igl_inline.h"
 
 #include <Eigen/Core>
 
-namespace igl {
-  /// Given 2 vectors centered on origin calculate the rotation matrix from first to the second
-
+namespace igl 
+{
+  // Given 2 vectors centered on origin calculate the rotation matrix from
+  // first to the second
+  //
   // Inputs:
-  //   v0, v1         the two #3 by 1 vectors
-  //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
+  //   v0  3D column vector
+  //   v1  3D column vector
   // Output:
-  //                  3 by 3 rotation matrix that takes v0 to v1
+  //   3 by 3 rotation matrix that takes v0 to v1
   //
   template <typename Scalar>
-  IGL_INLINE Eigen::Matrix<Scalar, 3, 3> rotation_matrix_from_directions(const Eigen::Matrix<Scalar, 3, 1> v0,
-                                                                     const Eigen::Matrix<Scalar, 3, 1> v1,
-                                                                     const bool normalized=true);
+  IGL_INLINE Eigen::Matrix<Scalar, 3, 3> rotation_matrix_from_directions(
+    const Eigen::Matrix<Scalar, 3, 1> v0,
+    const Eigen::Matrix<Scalar, 3, 1> v1);
 }
 
-
 #ifndef IGL_STATIC_LIBRARY
 #include "rotation_matrix_from_directions.cpp"
 #endif
-
-
-#endif /* defined(IGL_ROTATION_MATRIX_FROM_DIRECTIONS) */
+#endif

include/igl/sort.cpp

@@ -217,7 +217,7 @@ template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, 1, 3, 1, 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, 1, 1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
-template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
+template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<size_t,class std::allocator<size_t> > &);
 
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
@@ -237,8 +237,6 @@ template void igl::sort<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<do
 template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif
+template void igl::sort<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 
-#ifdef _WIN32
-template void igl::sort<int>(std::vector<int,std::allocator<int> > const &,bool,std::vector<int,std::allocator<int> > &,std::vector<unsigned long long int,class std::allocator<unsigned long long int> > &);
 #endif

include/igl/unproject_in_mesh.h

@@ -49,8 +49,9 @@ namespace igl
   //    model      model matrix
   //    proj       projection matrix
   //    viewport   vieweport vector
-  //    shoot_ray  function handle that outputs hits of a given ray against a
-  //      mesh (embedded in function handles as captured variable/data)
+  //    shoot_ray  function handle that outputs first hit of a given ray
+  //      against a mesh (embedded in function handles as captured
+  //      variable/data)
   // Outputs:
   //    obj        3d unprojected mouse point in mesh
   //    hits       vector of hits

include/igl/unproject_onto_mesh.cpp

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson
+//
+// 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 "unproject_onto_mesh.h"
+#include "unproject.h"
+#include "unproject_ray.h"
+#include "ray_mesh_intersect.h"
+#include <vector>
+
+template < typename DerivedV, typename DerivedF, typename Derivedbc>
+IGL_INLINE bool igl::unproject_onto_mesh(
+  const Eigen::Vector2f& pos,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  int & fid,
+  Eigen::PlainObjectBase<Derivedbc> & bc)
+{
+  using namespace std;
+  using namespace Eigen;
+  const auto & shoot_ray = [&V,&F](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    igl::Hit & hit)->bool
+  {
+    std::vector<igl::Hit> hits;
+    if(!ray_mesh_intersect(s,dir,V,F,hits))
+    {
+      return false;
+    }
+    hit = hits[0];
+    return true;
+  };
+  return unproject_onto_mesh(pos,model,proj,viewport,shoot_ray,fid,bc);
+}
+
+template <typename Derivedbc>
+IGL_INLINE bool igl::unproject_onto_mesh(
+  const Eigen::Vector2f& pos,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const std::function<
+    bool(
+      const Eigen::Vector3f&,
+      const Eigen::Vector3f&,
+      igl::Hit &)
+      > & shoot_ray,
+  int & fid,
+  Eigen::PlainObjectBase<Derivedbc> & bc)
+{
+  using namespace std;
+  using namespace Eigen;
+  Vector3f s,dir;
+  unproject_ray(pos,model,proj,viewport,s,dir);
+  Hit hit;
+  if(!shoot_ray(s,dir,hit))
+  {
+    return false;
+  }
+  bc.resize(3);
+  bc << 1.0-hit.u-hit.v, hit.u, hit.v;
+  fid = hit.id;
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::unproject_onto_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
+#endif
+

include/igl/unproject_onto_mesh.h

@@ -0,0 +1,74 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_UNPROJECT_ONTO_MESH
+#define IGL_UNPROJECT_ONTO_MESH
+#include "igl_inline.h"
+#include "Hit.h"
+#include <Eigen/Core>
+#include <functional>
+
+namespace igl
+{
+  // Unproject a screen location (using current opengl viewport, projection, and
+  // model view) to a 3D position _onto_ a given mesh, if the ray through the
+  // given screen location (x,y) _hits_ the mesh.
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    V   #V by 3 list of mesh vertex positions
+  //    F   #F by 3 list of mesh triangle indices into V
+  // Outputs:
+  //    fid  id of the first face hit
+  //    bc  barycentric coordinates of hit
+  // Returns true if there's a hit
+  template < typename DerivedV, typename DerivedF, typename Derivedbc>
+  IGL_INLINE bool unproject_onto_mesh(
+    const Eigen::Vector2f& pos,
+    const Eigen::Matrix4f& model,
+    const Eigen::Matrix4f& proj,
+    const Eigen::Vector4f& viewport,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    int & fid,
+    Eigen::PlainObjectBase<Derivedbc> & bc);
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  // Outputs:
+  //    fid  id of the first face hit
+  //    bc  barycentric coordinates of hit
+  // Returns true if there's a hit
+  template <typename Derivedbc>
+  IGL_INLINE bool unproject_onto_mesh(
+    const Eigen::Vector2f& pos,
+    const Eigen::Matrix4f& model,
+    const Eigen::Matrix4f& proj,
+    const Eigen::Vector4f& viewport,
+    const std::function<
+      bool(
+        const Eigen::Vector3f&,
+        const Eigen::Vector3f&,
+        igl::Hit  &)
+        > & shoot_ray,
+    int & fid,
+    Eigen::PlainObjectBase<Derivedbc> & bc);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "unproject_onto_mesh.cpp"
+#endif
+#endif
+
+

include/igl/viewer/Viewer.cpp

@@ -306,6 +306,22 @@ namespace viewer
     callback_mouse_scroll_data  = nullptr;
     callback_key_down_data      = nullptr;
     callback_key_up_data        = nullptr;
+
+#ifndef IGL_VIEWER_VIEWER_QUIET
+    const std::string usage(R"(igl::viewer::Viewer usage:
+  [drag]  Rotate scene
+  A,a     Toggle animation (tight draw loop)
+  I,i     Toggle invert normals
+  L,l     Toggle wireframe
+  O,o     Toggle orthographic/perspective projection
+  T,t     Toggle filled faces
+  Z       Snap to canonical view
+  [,]     Toggle between rotation control types (e.g. trackball, two-axis
+          valuator with fixed up)
+
+)");
+    std::cout<<usage;
+#endif
   }
 
   IGL_INLINE void Viewer::init_plugins()
@@ -456,18 +472,13 @@ namespace viewer
         return true;
 
     for (unsigned int i = 0; i<plugins.size(); ++i)
+    {
       if (plugins[i]->key_pressed(unicode_key, modifiers))
+      {
         return true;
-    const std::string usage(R"(igl::viewer::Viewer:
-A,a  Toggle animation (tight draw loop)
-I,i  Toggle invert normals
-L,l  Toggle wireframe
-O,o  Toggle orthographic/perspective projection
-T,t  Toggle filled faces
-Z    Snap to canonical view
-[,]  Toggle between rotation control types (e.g. trackball, two-axis valuator
-     with fixed up)
-)");
+      }
+    }
+
     switch(unicode_key)
     {
       case 'A':

include/igl/writePLY.cpp

@@ -132,7 +132,6 @@ IGL_INLINE bool igl::writePLY(
   }
 
   ply_close(ply);
-  fclose(fp);
   for(size_t i = 0;i<(size_t)F.rows();i++)
   {
     delete[] flist[i].verts;

tutorial/603_MEX/readOBJ_mex.cpp

@@ -36,9 +36,9 @@ void mexFunction(
   switch(nlhs)
   {
     case 2:
-      igl::prepare_lhs_index(F,plhs+1);
+      igl::matlab::prepare_lhs_index(F,plhs+1);
     case 1:
-      igl::prepare_lhs_double(V,plhs);
+      igl::matlab::prepare_lhs_double(V,plhs);
     default: break;
   }
 

tutorial/607_Picking/main.cpp

@@ -1,74 +1,50 @@
+#include "tutorial_shared_path.h"
 #include <igl/readOFF.h>
+#include <igl/unproject_onto_mesh.h>
 #include <igl/viewer/Viewer.h>
-#include <igl/embree/EmbreeIntersector.h>
-#include <igl/embree/unproject_onto_mesh.h>
-
-#include <algorithm>
-
-#include "tutorial_shared_path.h"
-
-// Mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-// Per-vertex color
-Eigen::MatrixXd C;
-
-igl::embree::EmbreeIntersector* ei;
-
-std::vector<int> picked_vertices;
-
-bool mouse_down(igl::viewer::Viewer& viewer, int button, int modifier)
-{
-
-  using namespace Eigen;
-  using namespace std;
-
-  int vid, fid;
-
-  // Cast a ray in the view direction starting from the mouse position
-  double x = viewer.current_mouse_x;
-  double y = viewer.core.viewport(3) - viewer.current_mouse_y;
-  bool hit = igl::embree::unproject_onto_mesh(Vector2f(x,y),
-                                F,
-                                viewer.core.view * viewer.core.model,
-                                viewer.core.proj,
-                                viewer.core.viewport,
-                                *ei,
-                                fid,
-                                vid);
-
-  // If the ray hits a face, assign a red color to the closest vertex
-  if (hit)
-  {
-    cerr << "Picked face(vertex): " << fid << " (" << vid << ")" << endl;
-    C.row(vid) << 1,0,0;
-    viewer.data.set_colors(C);
-    return true;
-  }
-
-  return false;
-}
+#include <iostream>
 
 int main(int argc, char *argv[])
 {
-
-  using namespace Eigen;
-  using namespace std;
+  // Mesh with per-face color
+  Eigen::MatrixXd V, C;
+  Eigen::MatrixXi F;
   // Load a mesh in OFF format
   igl::readOFF(TUTORIAL_SHARED_PATH "/fertility.off", V, F);
-
-  // Create a BVH for raycasting
-  ei = new igl::embree::EmbreeIntersector();
-  ei->init(V.cast<float>(),F);
-
-  // Initialize the colors to white for all vertices
-  C = MatrixXd::Constant(V.rows(),3,1);
-
-  // Show mesh
+  // Initialize white
+  C = Eigen::MatrixXd::Constant(F.rows(),3,1);
   igl::viewer::Viewer viewer;
+  viewer.callback_mouse_down = 
+    [&V,&F,&C](igl::viewer::Viewer& viewer, int, int)->bool
+  {
+    int fid;
+    Eigen::Vector3f bc;
+    // Cast a ray in the view direction starting from the mouse position
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    if(igl::unproject_onto_mesh(
+      Eigen::Vector2f(x,y),
+      viewer.core.view * viewer.core.model,
+      viewer.core.proj,
+      viewer.core.viewport,
+      V,
+      F,
+      fid,
+      bc))
+    {
+      // paint hit red
+      C.row(fid)<<1,0,0;
+      viewer.data.set_colors(C);
+      return true;
+    }
+    return false;
+  };
+  std::cout<<R"(Usage:
+  [click]  Pick face on shape
+
+)";
+  // Show mesh
   viewer.data.set_mesh(V, F);
-  viewer.callback_mouse_down = &mouse_down;
   viewer.data.set_colors(C);
   viewer.core.show_lines = false;
   viewer.launch();