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

Former-commit-id: 430e7d064217b89a79712821f5a1539460697ef1

include/igl/boolean/BinaryWindingNumberOperations.h

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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_BINARY_WINDING_NUMBER_OPERATIONS
+#define IGL_BINARY_WINDING_NUMBER_OPERATIONS
+
+#include <stdexcept>
+#include <igl/igl_inline.h>
+#include "MeshBooleanType.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+    namespace boolean
+    {
+        template <igl::boolean::MeshBooleanType Op>
+        class BinaryWindingNumberOperations {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        throw (std::runtime_error("not implemented!"));
+                    }
+        };
+
+        template <>
+        class BinaryWindingNumberOperations<igl::boolean::MESH_BOOLEAN_TYPE_UNION> {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        return win_nums(0) > 0 || win_nums(1) > 0;
+                    }
+        };
+
+        template <>
+        class BinaryWindingNumberOperations<igl::boolean::MESH_BOOLEAN_TYPE_INTERSECT> {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        return win_nums(0) > 0 && win_nums(1) > 0;
+                    }
+        };
+
+        template <>
+        class BinaryWindingNumberOperations<igl::boolean::MESH_BOOLEAN_TYPE_MINUS> {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        return win_nums(0) > 0 && win_nums(1) <= 0;
+                    }
+        };
+
+        template <>
+        class BinaryWindingNumberOperations<igl::boolean::MESH_BOOLEAN_TYPE_XOR> {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        return (win_nums(0) > 0 && win_nums(1) <= 0) ||
+                            (win_nums(0) <= 0 && win_nums(1) > 0);
+                    }
+        };
+
+        template <>
+        class BinaryWindingNumberOperations<igl::boolean::MESH_BOOLEAN_TYPE_RESOLVE> {
+            public:
+                template<typename DerivedW>
+                    typename DerivedW::Scalar operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        return true;
+                    }
+        };
+
+        typedef BinaryWindingNumberOperations<MESH_BOOLEAN_TYPE_UNION> BinaryUnion;
+        typedef BinaryWindingNumberOperations<MESH_BOOLEAN_TYPE_INTERSECT> BinaryIntersect;
+        typedef BinaryWindingNumberOperations<MESH_BOOLEAN_TYPE_MINUS> BinaryMinus;
+        typedef BinaryWindingNumberOperations<MESH_BOOLEAN_TYPE_XOR> BinaryXor;
+        typedef BinaryWindingNumberOperations<MESH_BOOLEAN_TYPE_RESOLVE> BinaryResolve;
+
+        enum KeeperType {
+            KEEP_INSIDE,
+            KEEP_ALL
+        };
+
+        template<KeeperType T>
+        class WindingNumberFilter {
+            public:
+                template<typename DerivedW>
+                    short operator()(
+                            const Eigen::PlainObjectBase<DerivedW>& win_nums) const {
+                        throw std::runtime_error("Not implemented");
+                    }
+        };
+
+        template<>
+        class WindingNumberFilter<KEEP_INSIDE> {
+            public:
+                template<typename T>
+                    short operator()(T out_w, T in_w) const {
+                        if (in_w > 0 && out_w <= 0) return 1;
+                        else if (in_w <= 0 && out_w > 0) return -1;
+                        else return 0;
+                    }
+        };
+
+        template<>
+        class WindingNumberFilter<KEEP_ALL> {
+            public:
+                template<typename T>
+                    short operator()(T out_w, T in_w) const {
+                        return 1;
+                    }
+        };
+
+        typedef WindingNumberFilter<KEEP_INSIDE> KeepInside;
+        typedef WindingNumberFilter<KEEP_ALL> KeepAll;
+    }
+}
+
+#endif

include/igl/boolean/CSGTree.h

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

include/igl/boolean/mesh_boolean.cpp

@@ -1,404 +1,369 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingnan Zhou <qnzhou@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 "mesh_boolean.h"
 #include <igl/cgal/assign_scalar.h>
-#include <igl/per_face_normals.h>
-#include <igl/boundary_facets.h>
-#include <igl/exterior_edges.h>
-#include <igl/cgal/peel_outer_hull_layers.h>
+#include <igl/cgal/propagate_winding_numbers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
-#include <igl/mod.h>
 #include <igl/unique_simplices.h>
+
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <iostream>
 
-//#define IGL_MESH_BOOLEAN_DEBUG
-template <
-  typename DerivedVA,
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename DerivedVC,
-  typename DerivedFC,
-  typename DerivedJ,
-  typename DerivedI>
-IGL_INLINE void igl::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  const std::function<void(
-    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
-    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-    Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
-    Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-    Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&,
-    Eigen::Matrix<typename DerivedI::Scalar, Eigen::Dynamic,1>&)> &
-    resolve_fun,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC,
-  Eigen::PlainObjectBase<DerivedJ > & J,
-  Eigen::PlainObjectBase<DerivedI > & I)
-{
-  using namespace Eigen;
-  using namespace std;
-  using namespace igl;
-  using namespace igl::cgal;
-  MeshBooleanType eff_type = type;
-  // Concatenate A and B into a single mesh
-  typedef CGAL::Epeck Kernel;
-  typedef Kernel::FT ExactScalar;
-  typedef typename DerivedVC::Scalar Scalar;
-  typedef typename DerivedFC::Scalar Index;
-  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
-  typedef Matrix<ExactScalar,Dynamic,3> MatrixX3ES;
-  typedef Matrix<Index,Dynamic,3> MatrixX3I;
-  typedef Matrix<Index,Dynamic,2> MatrixX2I;
-  typedef Matrix<typename DerivedI::Scalar,Dynamic,1> VectorXI;
-  typedef Matrix<typename DerivedJ::Scalar,Dynamic,1> VectorXJ;
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"mesh boolean..."<<endl;
-#endif
-  MatrixX3S V(VA.rows()+VB.rows(),3);
-  MatrixX3I F(FA.rows()+FB.rows(),3);
-  V.block(0,0,VA.rows(),VA.cols()) = VA;
-  V.block(VA.rows(),0,VB.rows(),VB.cols()) = VB;
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"prepare selfintersect input..."<<endl;
-#endif
-  switch(type)
-  {
-    // Minus is implemented by flipping B and computing union
-    case MESH_BOOLEAN_TYPE_MINUS:
-      F.block(0,0,FA.rows(),FA.cols()) = FA.rowwise().reverse();
-      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
-      //F.block(0,0,FA.rows(),3) = FA;
-      //F.block(FA.rows(),0,FB.rows(),3) =
-      //  FB.rowwise().reverse().array()+VA.rows();
-      eff_type = MESH_BOOLEAN_TYPE_INTERSECT;
-      break;
-    default:
-      F.block(0,0,FA.rows(),FA.cols()) = FA;
-      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
-      break;
-  }
+#include "BinaryWindingNumberOperations.h"
 
-  // Resolve intersections (assumes A and B are solid)
-  const auto & libigl_resolve = [](
-    const MatrixX3S & V,
-    const MatrixX3I & F,
-    MatrixX3ES & CV,
-    MatrixX3I & CF,
-    VectorXJ & J,
-    VectorXI & I)
-  {
-    MatrixX3ES SV;
-    MatrixX3I SF;
-    MatrixX2I SIF;
-    igl::cgal::RemeshSelfIntersectionsParam params;
-    remesh_self_intersections(V,F,params,SV,SF,SIF,J,I);
-    for_each(SF.data(),SF.data()+SF.size(),
-      [&I](typename MatrixX3I::Scalar & a){a=I(a);});
-    {
-      Eigen::Matrix<typename MatrixX3S::Index,Dynamic,1> UIM;
-      remove_unreferenced(SV,SF,CV,CF,UIM);
-      for_each(I.data(),I.data()+I.size(),
-        [&UIM](typename VectorXI::Scalar & a){a=UIM(a);});
-    }
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-    cout<<"#F:  "<<F.rows()<<endl;
-    cout<<"#CF: "<<CF.rows()<<endl;
-#endif
-  };
+namespace igl {
+    namespace boolean {
+        namespace mesh_boolean_helper {
+            typedef CGAL::Epeck Kernel;
+            typedef Kernel::FT ExactScalar;
 
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"resolve..."<<endl;
-#endif
-  MatrixX3S CV;
-  MatrixX3ES EV;
-  MatrixX3I CF;
-  VectorXJ CJ;
-  Eigen::Matrix<typename DerivedI::Scalar,Dynamic, 1> CI;
-  if(resolve_fun)
-  {
-    resolve_fun(V,F,CV,CF,CJ,CI);
-  }else
-  {
-    libigl_resolve(V,F,EV,CF,CJ,CI);
-    CV.resize(EV.rows(), EV.cols());
-    // Just use f'ing for loops. What if EV and CV don't use the same ordering?
-    for(int i=0;i<EV.rows();i++)
-    {
-      for(int j=0;j<EV.cols();j++)
-      {
-        assign_scalar(EV(i,j),CV(i,j));
-      }
-    }
-  }
+            template<
+                typename DerivedV,
+                typename DerivedF,
+                typename DerivedVo,
+                typename DerivedFo,
+                typename DerivedJ>
+            void igl_resolve(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    Eigen::PlainObjectBase<DerivedVo>& Vo,
+                    Eigen::PlainObjectBase<DerivedFo>& Fo,
+                    Eigen::PlainObjectBase<DerivedJ>& J) {
+                Eigen::VectorXi I;
+                igl::cgal::RemeshSelfIntersectionsParam params;
 
-  if(type == MESH_BOOLEAN_TYPE_RESOLVE)
-  {
-    FC = CF;
-    VC = CV;
-    J = CJ;
-    I = CI;
-    return;
-  }
+                DerivedVo Vr;
+                DerivedFo Fr;
+                Eigen::MatrixXi IF;
+                igl::cgal::remesh_self_intersections(V, F, params, Vr, Fr, IF, J, I);
+                assert(I.size() == Vr.rows());
 
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"peel..."<<endl;
-#endif
-  Matrix<bool,Dynamic,1> from_A(CF.rows());
-  // peel layers keeping track of odd and even flips
-  VectorXi iter;
-  Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(EV,CF,iter,flip);
-  //Array<bool,Dynamic,1> even = igl::mod(I,2).array()==0;
-  const auto even = [&](const Index & f)->bool
-  {
-    return (iter(f)%2)==0;
-  };
+                // Merge coinciding vertices into non-manifold vertices.
+                std::for_each(Fr.data(), Fr.data()+Fr.size(),
+                        [&I](typename DerivedF::Scalar& a) { a=I[a]; });
 
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"categorize..."<<endl;
-#endif
-  const Index m = CF.rows();
-  // Faces of output vG[i] = j means ith face of output should be jth face in F
-  std::vector<Index> vG;
-  // Whether faces of output should be flipped, Gflip[i] = true means ith face
-  // of output should be F.row(vG[i]).reverse() rather than F.row(vG[i])
-  std::vector<bool> Gflip;
-  for(Index f = 0;f<m;f++)
-  {
-    switch(eff_type)
-    {
-      case MESH_BOOLEAN_TYPE_XOR:
-      case MESH_BOOLEAN_TYPE_UNION:
-        if((even(f)&&!flip(f))||(!even(f)&&flip(f)))
-        {
-          vG.push_back(f);
-          Gflip.push_back(false);
-        }else if(eff_type == MESH_BOOLEAN_TYPE_XOR)
-        {
-          vG.push_back(f);
-          Gflip.push_back(true);
-        }
-        break;
-      case MESH_BOOLEAN_TYPE_INTERSECT:
-        if((!even(f) && !flip(f)) || (even(f) && flip(f)))
-        {
-          vG.push_back(f);
-          Gflip.push_back(type == MESH_BOOLEAN_TYPE_MINUS);
+                // Remove unreferenced vertices.
+                Eigen::VectorXi UIM;
+                igl::remove_unreferenced(Vr, Fr, Vo, Fo, UIM);
+            }
+
+            // Combine mesh A with mesh B and resolve all intersections.
+            template<
+                typename DerivedVA,
+                typename DerivedVB,
+                typename DerivedFA,
+                typename DerivedFB,
+                typename ResolveFunc,
+                typename DerivedVC,
+                typename DerivedFC,
+                typename DerivedJ>
+            void resolve_intersections(
+                    const Eigen::PlainObjectBase<DerivedVA>& VA,
+                    const Eigen::PlainObjectBase<DerivedFA>& FA,
+                    const Eigen::PlainObjectBase<DerivedVB>& VB,
+                    const Eigen::PlainObjectBase<DerivedFB>& FB,
+                    const ResolveFunc& resolve_func,
+                    Eigen::PlainObjectBase<DerivedVC>& VC,
+                    Eigen::PlainObjectBase<DerivedFC>& FC,
+                    Eigen::PlainObjectBase<DerivedJ>& J) {
+                DerivedVA V(VA.rows()+VB.rows(),3);
+                DerivedFA F(FA.rows()+FB.rows(),3);
+                V << VA, VB;
+                F << FA, FB.array() + VA.rows();
+                resolve_func(V, F, VC, FC, J);
+            }
+
+            template<
+                typename DerivedF1,
+                typename DerivedJ1,
+                typename DerivedF2,
+                typename DerivedJ2 >
+            void resolve_duplicated_faces(
+                    const Eigen::PlainObjectBase<DerivedF1>& F1,
+                    const Eigen::PlainObjectBase<DerivedJ1>& J1,
+                    Eigen::PlainObjectBase<DerivedF2>& F2,
+                    Eigen::PlainObjectBase<DerivedJ2>& J2) {
+                typedef typename DerivedF1::Scalar Index;
+                Eigen::VectorXi IA,IC;
+                DerivedF1 uF;
+                igl::unique_simplices(F1,uF,IA,IC);
+
+                const size_t num_faces = F1.rows();
+                const size_t num_unique_faces = uF.rows();
+                assert(IA.rows() == num_unique_faces);
+                // faces ontop of each unique face
+                std::vector<std::vector<int> > uF2F(num_unique_faces);
+                // signed counts
+                Eigen::VectorXi counts = Eigen::VectorXi::Zero(num_unique_faces);
+                Eigen::VectorXi ucounts = Eigen::VectorXi::Zero(num_unique_faces);
+                // loop over all faces
+                for (size_t i=0; i<num_faces; i++) {
+                    const size_t ui = IC(i);
+                    const bool consistent = 
+                        (F1(i,0) == uF(ui, 0) &&
+                         F1(i,1) == uF(ui, 1) &&
+                         F1(i,2) == uF(ui, 2)) ||
+                        (F1(i,0) == uF(ui, 1) &&
+                         F1(i,1) == uF(ui, 2) &&
+                         F1(i,2) == uF(ui, 0)) ||
+                        (F1(i,0) == uF(ui, 2) &&
+                         F1(i,1) == uF(ui, 0) &&
+                         F1(i,2) == uF(ui, 1));
+                    uF2F[ui].push_back(int(i+1) * (consistent?1:-1));
+                    counts(ui) += consistent ? 1:-1;
+                    ucounts(ui)++;
+                }
+
+                std::vector<size_t> kept_faces;
+                for (size_t i=0; i<num_unique_faces; i++) {
+                    if (ucounts[i] == 1) {
+                        kept_faces.push_back(abs(uF2F[i][0])-1);
+                        continue;
+                    }
+                    if (counts[i] == 1) {
+                        bool found = false;
+                        for (auto fid : uF2F[i]) {
+                            if (fid > 0) {
+                                kept_faces.push_back(abs(fid)-1);
+                                found = true;
+                                break;
+                            }
+                        }
+                        assert(found);
+                    } else if (counts[i] == -1) {
+                        bool found = false;
+                        for (auto fid : uF2F[i]) {
+                            if (fid < 0) {
+                                kept_faces.push_back(abs(fid)-1);
+                                found = true;
+                                break;
+                            }
+                        }
+                        assert(found);
+                    } else {
+                        assert(counts[i] == 0);
+                    }
+                }
+
+                const size_t num_kept = kept_faces.size();
+                F2.resize(num_kept, 3);
+                J2.resize(num_kept, 1);
+                for (size_t i=0; i<num_kept; i++) {
+                    F2.row(i) = F1.row(kept_faces[i]);
+                    J2.row(i) = J1.row(kept_faces[i]);
+                }
+            }
         }
-        break;
-      default:
-        assert(false && "Unknown type");
-        return;
     }
-  }
-  const Index gm = vG.size();
-  MatrixX3I G(gm,3);
-  VectorXi GJ(gm,1);
-  for(Index g = 0;g<gm;g++)
-  {
-    G.row(g) = Gflip[g] ? CF.row(vG[g]).reverse().eval() : CF.row(vG[g]);
-    GJ(g) = CJ(vG[g]);
-  }
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  {
-    MatrixXi O;
-    boundary_facets(FC,O);
-    cout<<"# boundary: "<<O.rows()<<endl;
-  }
-  cout<<"# exterior: "<<exterior_edges(FC).rows()<<endl;
-#endif
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  cout<<"clean..."<<endl;
-#endif
-  // Deal with duplicate faces
-  {
-    VectorXi IA,IC;
-    MatrixX3I uG;
-    unique_simplices(G,uG,IA,IC);
-    assert(IA.rows() == uG.rows());
-    // faces ontop of each unique face
-    vector<vector<Index> > uG2G(uG.rows());
-    // signed counts
-    VectorXi counts = VectorXi::Zero(uG.rows());
-    VectorXi ucounts = VectorXi::Zero(uG.rows());
-    // loop over all faces
-    for(Index g = 0;g<gm;g++)
-    {
-      const int ug = IC(g);
-      assert((size_t) ug < uG2G.size());
-      uG2G[ug].push_back(g);
-      // is uG(g,:) just a rotated version of G(g,:) ?
-      const bool consistent =
-        (G(g,0) == uG(ug,0) && G(g,1) == uG(ug,1) && G(g,2) == uG(ug,2)) ||
-        (G(g,0) == uG(ug,1) && G(g,1) == uG(ug,2) && G(g,2) == uG(ug,0)) ||
-        (G(g,0) == uG(ug,2) && G(g,1) == uG(ug,0) && G(g,2) == uG(ug,1));
-      counts(ug) += consistent ? 1 : -1;
-      ucounts(ug)++;
+}
+
+template <
+typename DerivedVA,
+typename DerivedFA,
+typename DerivedVB,
+typename DerivedFB,
+typename WindingNumberOp,
+typename KeepFunc,
+typename ResolveFunc,
+typename DerivedVC,
+typename DerivedFC,
+typename DerivedJ>
+IGL_INLINE void igl::boolean::per_face_winding_number_binary_operation(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        const WindingNumberOp& wind_num_op,
+        const KeepFunc& keep,
+        const ResolveFunc& resolve_fun,
+        Eigen::PlainObjectBase<DerivedVC > & VC,
+        Eigen::PlainObjectBase<DerivedFC > & FC,
+        Eigen::PlainObjectBase<DerivedJ > & J) {
+    using namespace igl::boolean::mesh_boolean_helper;
+
+    typedef typename DerivedVC::Scalar Scalar;
+    typedef typename DerivedFC::Scalar Index;
+    typedef Eigen::Matrix<Scalar,Eigen::Dynamic,3> MatrixX3S;
+    typedef Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic> MatrixXI;
+    typedef Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1> VectorXJ;
+
+    // Generate combined mesh.
+    typedef Eigen::Matrix<
+        ExactScalar,
+        Eigen::Dynamic,
+        Eigen::Dynamic,
+        DerivedVC::IsRowMajor> MatrixXES;
+    MatrixXES V;
+    DerivedFC F;
+    VectorXJ  CJ;
+    resolve_intersections(VA, FA, VB, FB, resolve_fun, V, F, CJ);
+
+    // Compute winding numbers on each side of each facet.
+    const size_t num_faces = F.rows();
+    Eigen::MatrixXi W;
+    Eigen::VectorXi labels(num_faces);
+    std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
+            [&](int i) { return i<FA.rows() ? 0:1; });
+    igl::cgal::propagate_winding_numbers(V, F, labels, W);
+    assert(W.rows() == num_faces);
+    if (W.cols() == 2) {
+        assert(FB.rows() == 0);
+        Eigen::MatrixXi W_tmp(num_faces, 4);
+        W_tmp << W, Eigen::MatrixXi::Zero(num_faces, 2);
+        W = W_tmp;
+    } else {
+        assert(W.cols() == 4);
     }
-    MatrixX3I oldG = G;
-    // Faces of output vG[i] = j means ith face of output should be jth face in
-    // oldG
-    vG.clear();
-    for(size_t ug = 0;ug < uG2G.size();ug++)
-    {
-      // if signed occurrences is zero or ±two then keep none
-      // else if signed occurrences is ±one then keep just one facet
-      switch(abs(counts(ug)))
-      {
-        case 1:
-          assert(uG2G[ug].size() > 0);
-          vG.push_back(uG2G[ug][0]);
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-          if(abs(ucounts(ug)) != 1)
-          {
-            cout<<"count,ucount of "<<counts(ug)<<","<<ucounts(ug)<<endl;
-          }
-#endif
-          break;
-        case 0:
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-          cout<<"Skipping "<<uG2G[ug].size()<<" facets..."<<endl;
-          if(abs(ucounts(ug)) != 0)
-          {
-            cout<<"count,ucount of "<<counts(ug)<<","<<ucounts(ug)<<endl;
-          }
-#endif
-          break;
-        default:
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-          cout<<"Didn't expect to be here."<<endl;
-#endif
-          assert(false && "Shouldn't count be -1/0/1 ?");
-      }
+
+    // Compute resulting winding number.
+    Eigen::MatrixXi Wr(num_faces, 2);
+    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);
+        Wr(i,0) = wind_num_op(w_out);
+        Wr(i,1) = wind_num_op(w_in);
     }
-    G.resize(vG.size(),3);
-    J.resize(vG.size());
-    for(size_t g = 0;g<vG.size();g++)
-    {
-      G.row(g) = oldG.row(vG[g]);
-      J(g) = GJ(vG[g]);
+
+    // Extract boundary separating inside from outside.
+    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++) {
+        auto should_keep = keep(Wr(i,0), Wr(i,1));
+        if (should_keep > 0) {
+            selected.push_back(index_to_signed_index(i, true));
+        } else if (should_keep < 0) {
+            selected.push_back(index_to_signed_index(i, false));
+        }
+    }
+
+    const size_t num_selected = selected.size();
+    DerivedFC kept_faces(num_selected, 3);
+    DerivedJ  kept_face_indices;
+    kept_face_indices.resize(num_selected, 1);
+    for (size_t i=0; i<num_selected; i++) {
+        size_t idx = abs(selected[i]) - 1;
+        if (selected[i] > 0) {
+            kept_faces.row(i) = F.row(idx);
+        } else {
+            kept_faces.row(i) = F.row(idx).reverse();
+        }
+        kept_face_indices(i, 0) = CJ[idx];
     }
-  }
-  // remove unreferenced vertices
-  {
-    VectorXi newIM;
-    remove_unreferenced(CV,G,VC,FC,newIM);
-    I.resize(CI.rows(),CI.cols());
-    for(int i = 0;i<CI.rows();i++)
+
+
+    // Finally, remove duplicated faces and unreferenced vertices.
     {
-      I(i) = newIM(CI(i));
+        DerivedFC G;
+        DerivedJ J;
+        resolve_duplicated_faces(kept_faces, kept_face_indices, G, J);
+
+        MatrixX3S Vs(V.rows(), V.cols());
+        for (size_t i=0; i<V.rows(); i++) {
+            for (size_t j=0; j<V.cols(); j++) {
+                igl::cgal::assign_scalar(V(i,j), Vs(i,j));
+            }
+        }
+        Eigen::VectorXi newIM;
+        igl::remove_unreferenced(Vs,G,VC,FC,newIM);
     }
-  }
-  //cerr<<"warning not removing unref"<<endl;
-  //VC = CV;
-  //FC = G;
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-  {
-    MatrixXi O;
-    boundary_facets(FC,O);
-    cout<<"# boundary: "<<O.rows()<<endl;
-  }
-  cout<<"# exterior: "<<exterior_edges(FC).rows()<<endl;
-#endif
 }
 
-
 template <
-  typename DerivedVA,
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename DerivedVC,
-  typename DerivedFC,
-  typename DerivedJ,
-  typename DerivedI>
+typename DerivedVA,
+typename DerivedFA,
+typename DerivedVB,
+typename DerivedFB,
+typename DerivedVC,
+typename DerivedFC,
+typename DerivedJ>
 IGL_INLINE void igl::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC,
-  Eigen::PlainObjectBase<DerivedJ > & J,
-  Eigen::PlainObjectBase<DerivedI > & I)
-{
-  const std::function<void(
-    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&,
-          Eigen::Matrix<typename DerivedI::Scalar, Eigen::Dynamic,1>&) >
-    empty_fun;
-  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J,I);
+        const Eigen::PlainObjectBase<DerivedVA > & VA,
+        const Eigen::PlainObjectBase<DerivedFA > & FA,
+        const Eigen::PlainObjectBase<DerivedVB > & VB,
+        const Eigen::PlainObjectBase<DerivedFB > & FB,
+        const MeshBooleanType & type,
+        Eigen::PlainObjectBase<DerivedVC > & VC,
+        Eigen::PlainObjectBase<DerivedFC > & FC,
+        Eigen::PlainObjectBase<DerivedJ > & J) {
+    using namespace igl::boolean::mesh_boolean_helper;
+    typedef Eigen::Matrix<
+        ExactScalar,
+        Eigen::Dynamic,
+        Eigen::Dynamic,
+        DerivedVC::IsRowMajor> MatrixXES;
+    std::function<void(
+            const Eigen::PlainObjectBase<DerivedVA>&,
+            const Eigen::PlainObjectBase<DerivedFA>&,
+            Eigen::PlainObjectBase<MatrixXES>&,
+            Eigen::PlainObjectBase<DerivedFC>&,
+            Eigen::PlainObjectBase<DerivedJ>&)> resolve_func =
+        igl_resolve<DerivedVA, DerivedFA, MatrixXES, DerivedFC, DerivedJ>;
+
+    switch (type) {
+        case MESH_BOOLEAN_TYPE_UNION:
+            igl::boolean::per_face_winding_number_binary_operation(
+                    VA, FA, VB, FB, igl::boolean::BinaryUnion(),
+                    igl::boolean::KeepInside(), resolve_func, VC, FC, J);
+            break;
+        case MESH_BOOLEAN_TYPE_INTERSECT:
+            igl::boolean::per_face_winding_number_binary_operation(
+                    VA, FA, VB, FB, igl::boolean::BinaryIntersect(),
+                    igl::boolean::KeepInside(), resolve_func, VC, FC, J);
+            break;
+        case MESH_BOOLEAN_TYPE_MINUS:
+            igl::boolean::per_face_winding_number_binary_operation(
+                    VA, FA, VB, FB, igl::boolean::BinaryMinus(),
+                    igl::boolean::KeepInside(), resolve_func, VC, FC, J);
+            break;
+        case MESH_BOOLEAN_TYPE_XOR:
+            igl::boolean::per_face_winding_number_binary_operation(
+                    VA, FA, VB, FB, igl::boolean::BinaryXor(),
+                    igl::boolean::KeepInside(), resolve_func, VC, FC, J);
+            break;
+        case MESH_BOOLEAN_TYPE_RESOLVE:
+            //op = binary_resolve();
+            igl::boolean::per_face_winding_number_binary_operation(
+                    VA, FA, VB, FB, igl::boolean::BinaryResolve(),
+                    igl::boolean::KeepAll(), resolve_func, VC, FC, J);
+            break;
+        default:
+            throw std::runtime_error("Unsupported boolean type.");
+    }
 }
 
 template <
-  typename DerivedVA,
-  typename DerivedFA,
-  typename DerivedVB,
-  typename DerivedFB,
-  typename DerivedVC,
-  typename DerivedFC>
+typename DerivedVA,
+typename DerivedFA,
+typename DerivedVB,
+typename DerivedFB,
+typename DerivedVC,
+typename DerivedFC>
 IGL_INLINE void igl::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC)
-{
-  Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
-  typedef Eigen::Matrix<typename DerivedVC::Index, Eigen::Dynamic,1> VectorXI;
-  VectorXI I;
-  const std::function<void(
-    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&,
-          Eigen::Matrix<typename VectorXI::Scalar, Eigen::Dynamic,1>&)>
-    empty_fun;
-  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J,I);
+        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::boolean::mesh_boolean(VA,FA,VB,FB,type,VC,FC,J);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::boolean::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-// This is a hack to discuss. I'm not sure why this _doesn't_ create
-// duplicate symbols.
-#include <igl/remove_unreferenced.cpp>
-template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-#include <igl/cgal/peel_outer_hull_layers.cpp>
-template unsigned long
-igl::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-#include <igl/cgal/outer_hull.cpp>
-template void igl::cgal::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-#include <igl/slice.cpp>
-#include <igl/barycenter.cpp>
-template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
-#include <igl/mod.cpp>
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > igl::mod<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int);
-#include <igl/outer_element.cpp>
-template void igl::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-#include <igl/colon.cpp>
-template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
-
-
 #endif

include/igl/boolean/mesh_boolean.h

@@ -1,10 +1,12 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingnan Zhou <qnzhou@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_BOOLEAN_MESH_BOOLEAN_H
 #define IGL_BOOLEAN_MESH_BOOLEAN_H
 
@@ -17,54 +19,29 @@ namespace igl
 {
   namespace boolean
   {
-    //  MESH_BOOLEAN Compute boolean csg operations on "solid", consistently
-    //  oriented meshes.
-    // 
-    //  Inputs:
-    //    VA  #VA by 3 list of vertex positions of first mesh
-    //    FA  #FA by 3 list of triangle indices into VA
-    //    VB  #VB by 3 list of vertex positions of second mesh
-    //    FB  #FB by 3 list of triangle indices into VB
-    //    type  type of boolean operation
-    //    resolve_fun  function handle for computing resolve of a
-    //      self-intersections of a mesh and outputting the new mesh.
-    //  Outputs:
-    //    VC  #VC by 3 list of vertex positions of boolean result mesh
-    //    FC  #FC by 3 list of triangle indices into VC
-    //    J  #FC list of indices into [FA;FB] revealing "birth" facet
-    //    I  #VA+#VB list of indices into SV (SV=[VA;VB;SVA;SVB], where SVA and
-    //      SVB are the new vertices from resolving intersections) revealing
-    //      "birth" vertices.
-    //
-    //  See also: mesh_boolean_cork, intersect_other, remesh_self_intersections
-    //     
-    template <
-      typename DerivedVA,
-      typename DerivedFA,
-      typename DerivedVB,
-      typename DerivedFB,
-      typename DerivedVC,
-      typename DerivedFC,
-      typename DerivedJ,
-      typename DerivedI>
-    IGL_INLINE void mesh_boolean(
-      const Eigen::PlainObjectBase<DerivedVA > & VA,
-      const Eigen::PlainObjectBase<DerivedFA > & FA,
-      const Eigen::PlainObjectBase<DerivedVB > & VB,
-      const Eigen::PlainObjectBase<DerivedFB > & FB,
-      const MeshBooleanType & type,
-      const std::function<void(
-        const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-        const Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
-        Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-        Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
-        Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1>&,
-        Eigen::Matrix<typename DerivedI::Scalar,Eigen::Dynamic,1>&)> &
-        resolve_fun,
-      Eigen::PlainObjectBase<DerivedVC > & VC,
-      Eigen::PlainObjectBase<DerivedFC > & FC,
-      Eigen::PlainObjectBase<DerivedJ > & J,
-      Eigen::PlainObjectBase<DerivedI > & I);
+      template <
+          typename DerivedVA,
+          typename DerivedFA,
+          typename DerivedVB,
+          typename DerivedFB,
+          typename WindingNumberOp,
+          typename IsInsideFunc,
+          typename ResolveFunc,
+          typename DerivedVC,
+          typename DerivedFC,
+          typename DerivedJ>
+      IGL_INLINE void per_face_winding_number_binary_operation(
+              const Eigen::PlainObjectBase<DerivedVA> & VA,
+              const Eigen::PlainObjectBase<DerivedFA> & FA,
+              const Eigen::PlainObjectBase<DerivedVB> & VB,
+              const Eigen::PlainObjectBase<DerivedFB> & FB,
+              const WindingNumberOp& wind_num_op,
+              const IsInsideFunc& is_inside,
+              const ResolveFunc& resolve_fun,
+              Eigen::PlainObjectBase<DerivedVC > & VC,
+              Eigen::PlainObjectBase<DerivedFC > & FC,
+              Eigen::PlainObjectBase<DerivedJ > & J);
+
     template <
       typename DerivedVA,
       typename DerivedFA,
@@ -72,8 +49,7 @@ namespace igl
       typename DerivedFB,
       typename DerivedVC,
       typename DerivedFC,
-      typename DerivedJ,
-      typename DerivedI>
+      typename DerivedJ>
     IGL_INLINE void mesh_boolean(
       const Eigen::PlainObjectBase<DerivedVA > & VA,
       const Eigen::PlainObjectBase<DerivedFA > & FA,
@@ -82,8 +58,8 @@ namespace igl
       const MeshBooleanType & type,
       Eigen::PlainObjectBase<DerivedVC > & VC,
       Eigen::PlainObjectBase<DerivedFC > & FC,
-      Eigen::PlainObjectBase<DerivedJ > & J,
-      Eigen::PlainObjectBase<DerivedI > & I);
+      Eigen::PlainObjectBase<DerivedJ > & J);
+
     template <
       typename DerivedVA,
       typename DerivedFA,

include/igl/cgal/SelfIntersectMesh.h

@@ -225,6 +225,7 @@ namespace igl
 
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
+#include "remesh_intersections.h"
 
 #include <igl/REDRUM.h>
 #include <igl/get_seconds.h>

include/igl/cgal/closest_facet.cpp

@@ -0,0 +1,356 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "closest_facet.h"
+
+#include <vector>
+#include <stdexcept>
+
+#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"
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI,
+    typename DerivedP,
+    typename DerivedR,
+    typename DerivedS >
+IGL_INLINE void igl::cgal::closest_facet(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedR>& R,
+        Eigen::PlainObjectBase<DerivedS>& S) {
+    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;
+
+    if (F.rows() <= 0 || I.rows() <= 0) {
+        throw std::runtime_error(
+                "Closest facet cannot be computed on empty mesh.");
+    }
+
+    const size_t num_faces = I.rows();
+    std::vector<Triangle> triangles;
+    for (size_t i=0; i<num_faces; i++) {
+        const Eigen::Vector3i f = F.row(I(i, 0));
+        triangles.emplace_back(
+                Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+                Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+                Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+        if (triangles.back().is_degenerate()) {
+            throw std::runtime_error(
+                    "Input facet components contains degenerated triangles");
+        }
+    }
+    Tree tree(triangles.begin(), triangles.end());
+    tree.accelerate_distance_queries();
+
+    auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
+        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));
+        Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+        auto ori = CGAL::orientation(v0, v1, v2, p);
+        switch (ori) {
+            case CGAL::POSITIVE:
+                return true;
+            case CGAL::NEGATIVE:
+                return false;
+            case CGAL::COPLANAR:
+                throw std::runtime_error(
+                        "It seems input mesh contains self intersection");
+            default:
+                throw std::runtime_error("Unknown CGAL state.");
+        }
+        return false;
+    };
+
+    auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
+        const auto& f = F.row(fid);
+        if (f[0] == s && f[1] == d) return false;
+        else if (f[1] == s && f[2] == d) return false;
+        else if (f[2] == s && f[0] == d) return false;
+        else if (f[0] == d && f[1] == s) return true;
+        else if (f[1] == d && f[2] == s) return true;
+        else if (f[2] == d && f[0] == s) return true;
+        else {
+            throw std::runtime_error(
+                    "Cannot compute orientation due to incorrect connectivity");
+            return false;
+        }
+    };
+    auto index_to_signed_index = [&](size_t index, bool ori) -> int{
+        return (index+1) * (ori? 1:-1);
+    };
+    auto signed_index_to_index = [&](int signed_index) -> size_t {
+        return abs(signed_index) - 1;
+    };
+
+    enum ElementType { VERTEX, EDGE, FACE };
+    auto determine_element_type = [&](const Point_3& p, const size_t fid,
+            size_t& element_index) {
+        const auto& tri = triangles[fid];
+        const Point_3 p0 = tri[0];
+        const Point_3 p1 = tri[1];
+        const Point_3 p2 = tri[2];
+
+        if (p == p0) { element_index = 0; return VERTEX; }
+        if (p == p1) { element_index = 1; return VERTEX; }
+        if (p == p2) { element_index = 2; return VERTEX; }
+        if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+        if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+        if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+        element_index = 0;
+        return FACE;
+    };
+
+    auto process_edge_case = [&](
+            size_t query_idx,
+            const size_t s, const size_t d,
+            size_t preferred_facet,
+            bool& orientation) {
+
+        Point_3 mid_edge_point(
+                (V(s,0) + V(d,0)) * 0.5,
+                (V(s,1) + V(d,1)) * 0.5,
+                (V(s,2) + V(d,2)) * 0.5);
+        Point_3 query_point(
+                P(query_idx, 0),
+                P(query_idx, 1),
+                P(query_idx, 2));
+
+        std::vector<Tree::Primitive_id> intersected_faces;
+        tree.all_intersected_primitives(Segment_3(mid_edge_point, query_point),
+                std::back_inserter(intersected_faces));
+
+        const size_t num_intersected_faces = intersected_faces.size();
+        std::vector<size_t> intersected_face_indices(num_intersected_faces);
+        std::vector<int> intersected_face_signed_indices(num_intersected_faces);
+        std::transform(intersected_faces.begin(),
+                intersected_faces.end(),
+                intersected_face_indices.begin(),
+                [&](const Tree::Primitive_id& itr) -> size_t
+                { return I(itr-triangles.begin(), 0); });
+        std::transform(
+                intersected_face_indices.begin(),
+                intersected_face_indices.end(),
+                intersected_face_signed_indices.begin(),
+                [&](size_t index) {
+                    return index_to_signed_index(
+                        index, get_orientation(index, s,d));
+                });
+
+        assert(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.
+            const size_t fid = intersected_face_indices[0];
+            orientation = on_the_positive_side(fid, query_point);
+            return fid;
+        }
+
+        Eigen::VectorXi order;
+        DerivedP pivot = P.row(query_idx).eval();
+        igl::cgal::order_facets_around_edge(V, F, s, d,
+                intersected_face_signed_indices,
+                pivot, order);
+
+        // Although first and last are equivalent, make the choice based on
+        // preferred_facet.
+        const size_t first = order[0];
+        const size_t last = order[num_intersected_faces-1];
+        if (intersected_face_indices[first] == preferred_facet) {
+            orientation = intersected_face_signed_indices[first] < 0;
+            return intersected_face_indices[first];
+        } else if (intersected_face_indices[last] == preferred_facet) {
+            orientation = intersected_face_signed_indices[last] > 0;
+            return intersected_face_indices[last];
+        } else {
+            orientation = intersected_face_signed_indices[order[0]] < 0;
+            return intersected_face_indices[order[0]];
+        }
+    };
+
+    auto process_face_case = [&](
+            const size_t query_idx, const size_t fid, bool& orientation) {
+        const auto& f = F.row(I(fid, 0));
+        return process_edge_case(query_idx, f[0], f[1], I(fid, 0), orientation);
+    };
+
+    auto process_vertex_case = [&](const size_t query_idx, size_t s,
+            size_t preferred_facet, bool& orientation) {
+        Point_3 closest_point(V(s, 0), V(s, 1), V(s, 2));
+        Point_3 query_point(P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
+
+        std::vector<Tree::Primitive_id> intersected_faces;
+        tree.all_intersected_primitives(Segment_3(closest_point, query_point),
+                std::back_inserter(intersected_faces));
+
+        const size_t num_intersected_faces = intersected_faces.size();
+        std::vector<size_t> intersected_face_indices(num_intersected_faces);
+        std::transform(intersected_faces.begin(),
+                intersected_faces.end(),
+                intersected_face_indices.begin(),
+                [&](const Tree::Primitive_id& itr) -> size_t
+                { return I(itr-triangles.begin(), 0); });
+
+        std::set<size_t> adj_vertices_set;
+        for (auto fid : intersected_face_indices) {
+            const auto& f = F.row(fid);
+            if (f[0] != s) adj_vertices_set.insert(f[0]);
+            if (f[1] != s) adj_vertices_set.insert(f[1]);
+            if (f[2] != s) adj_vertices_set.insert(f[2]);
+        }
+        const size_t num_adj_vertices = adj_vertices_set.size();
+        std::vector<size_t> adj_vertices(num_adj_vertices);
+        std::copy(adj_vertices_set.begin(), adj_vertices_set.end(),
+                adj_vertices.begin());
+
+        std::vector<Point_3> adj_points;
+        for (size_t vid : adj_vertices) {
+            adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
+        }
+
+        // A plane is on the exterior if all adj_points lies on or to
+        // one side of the plane.
+        auto is_on_exterior = [&](const Plane_3& separator) {
+            size_t positive=0;
+            size_t negative=0;
+            size_t coplanar=0;
+            for (const auto& point : adj_points) {
+                switch(separator.oriented_side(point)) {
+                    case CGAL::ON_POSITIVE_SIDE:
+                        positive++;
+                        break;
+                    case CGAL::ON_NEGATIVE_SIDE:
+                        negative++;
+                        break;
+                    case CGAL::ON_ORIENTED_BOUNDARY:
+                        coplanar++;
+                        break;
+                    default:
+                        throw "Unknown plane-point orientation";
+                }
+            }
+            auto query_orientation = separator.oriented_side(query_point);
+            bool r = (positive == 0 && query_orientation == CGAL::POSITIVE)
+                || (negative == 0 && query_orientation == CGAL::NEGATIVE);
+            return r;
+        };
+
+        size_t d = std::numeric_limits<size_t>::max();
+        for (size_t i=0; i<num_adj_vertices; i++) {
+            const size_t vi = adj_vertices[i];
+            for (size_t j=i+1; j<num_adj_vertices; j++) {
+                Plane_3 separator(closest_point, adj_points[i], adj_points[j]);
+                if (separator.is_degenerate()) {
+                    throw std::runtime_error(
+                            "Input mesh contains degenerated faces");
+                }
+                if (is_on_exterior(separator)) {
+                    d = vi;
+                    assert(!CGAL::collinear(
+                                query_point, adj_points[i], closest_point));
+                    break;
+                }
+            }
+        }
+        assert(d != std::numeric_limits<size_t>::max());
+
+        return process_edge_case(query_idx, s, d, preferred_facet, orientation);
+    };
+
+    const size_t num_queries = P.rows();
+    R.resize(num_queries, 1);
+    S.resize(num_queries, 1);
+    for (size_t i=0; i<num_queries; i++) {
+        const Point_3 query(P(i,0), P(i,1), P(i,2));
+        auto projection = tree.closest_point_and_primitive(query);
+        const Point_3 closest_point = projection.first;
+        size_t fid = projection.second - triangles.begin();
+        bool fid_ori = false;
+
+        // Gether all facets sharing the closest point.
+        std::vector<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();
+        std::vector<size_t> intersected_face_indices(num_intersected_faces);
+        std::transform(intersected_faces.begin(),
+                intersected_faces.end(),
+                intersected_face_indices.begin(),
+                [&](const Tree::Primitive_id& itr) -> size_t
+                { return I(itr-triangles.begin(), 0); });
+
+        size_t element_index;
+        auto element_type = determine_element_type(closest_point, fid,
+                element_index);
+        switch(element_type) {
+            case VERTEX:
+                {
+                    const auto& f = F.row(I(fid, 0));
+                    const size_t s = f[element_index];
+                    fid = process_vertex_case(i, s, I(fid, 0), fid_ori);
+                }
+                break;
+            case EDGE:
+                {
+                    const auto& f = F.row(I(fid, 0));
+                    const size_t s = f[(element_index+1)%3];
+                    const size_t d = f[(element_index+2)%3];
+                    fid = process_edge_case(i, s, d, I(fid, 0), fid_ori);
+                }
+                break;
+            case FACE:
+                {
+                    fid = process_face_case(i, fid, fid_ori);
+                }
+                break;
+            default:
+                throw std::runtime_error("Unknown element type.");
+        }
+
+
+        R(i,0) = fid;
+        S(i,0) = fid_ori;
+    }
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedR,
+    typename DerivedS >
+IGL_INLINE void igl::cgal::closest_facet(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedR>& R,
+        Eigen::PlainObjectBase<DerivedS>& S) {
+    const size_t num_faces = F.rows();
+    Eigen::VectorXi I(num_faces);
+    I.setLinSpaced(num_faces, 0, num_faces-1);
+    igl::cgal::closest_facet(V, F, I, P, R, S);
+}

include/igl/cgal/closest_facet.h

@@ -0,0 +1,63 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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 CLOSEST_FACET_H
+#define CLOSEST_FACET_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+    namespace cgal {
+
+        // Determine the closest facet for each of the input points.
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //   I  #I list of triangle indices to consider.
+        //   P  #P by 3 array of query points.
+        //
+        // Outputs:
+        //   R  #P list of closest facet indices.
+        //   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 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,
+                Eigen::PlainObjectBase<DerivedR>& R,
+                Eigen::PlainObjectBase<DerivedS>& S);
+
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedP,
+            typename DerivedR,
+            typename DerivedS >
+        IGL_INLINE void closest_facet(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                const Eigen::PlainObjectBase<DerivedP>& P,
+                Eigen::PlainObjectBase<DerivedR>& R,
+                Eigen::PlainObjectBase<DerivedS>& S);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "closest_facet.cpp"
+#endif
+#endif

include/igl/cgal/extract_cells.cpp

@@ -0,0 +1,340 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "extract_cells.h"
+#include "../extract_manifold_patches.h"
+#include "../facet_components.h"
+#include "../triangle_triangle_adjacency.h"
+#include "../unique_edge_map.h"
+#include "closest_facet.h"
+#include "order_facets_around_edge.h"
+#include "outer_facet.h"
+
+#include <vector>
+#include <queue>
+
+template<
+typename DerivedV,
+typename DerivedF,
+typename DerivedP,
+typename DeriveduE,
+typename uE2EType,
+typename DerivedEMAP,
+typename DerivedC >
+IGL_INLINE size_t igl::cgal::extract_cells_single_component(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        Eigen::PlainObjectBase<DerivedC>& cells) {
+    typedef typename DerivedF::Scalar Index;
+    const size_t num_faces = F.rows();
+    auto edge_index_to_face_index = [&](size_t index) {
+        return index % num_faces;
+    };
+    auto is_consistent = [&](const size_t fid, const size_t s, const size_t d) {
+        if (F(fid, 0) == s && F(fid, 1) == d) return false;
+        if (F(fid, 1) == s && F(fid, 2) == d) return false;
+        if (F(fid, 2) == s && F(fid, 0) == d) return false;
+
+        if (F(fid, 0) == d && F(fid, 1) == s) return true;
+        if (F(fid, 1) == d && F(fid, 2) == s) return true;
+        if (F(fid, 2) == d && F(fid, 0) == s) return true;
+        throw "Invalid face!";
+        return false;
+    };
+
+    const size_t num_unique_edges = uE.rows();
+    const size_t num_patches = P.maxCoeff() + 1;
+
+    std::vector<std::vector<size_t> > patch_edge_adj(num_patches);
+    std::vector<Eigen::VectorXi> orders(num_unique_edges);
+    std::vector<std::vector<bool> > orientations(num_unique_edges);
+    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 (adj_faces.size() > 2) {
+            std::vector<int> signed_adj_faces;
+            for (auto ei : adj_faces) {
+                const size_t fid = edge_index_to_face_index(ei);
+                bool cons = is_consistent(fid, s, d);
+                signed_adj_faces.push_back((fid+1)*(cons ? 1:-1));
+            }
+            auto& order = orders[i];
+            igl::cgal::order_facets_around_edge(
+                    V, F, s, d, signed_adj_faces, order);
+            auto& orientation = orientations[i];
+            orientation.resize(order.size());
+            std::transform(order.data(), order.data() + order.size(),
+                    orientation.begin(), [&](int index) { return
+                    signed_adj_faces[index] > 0; });
+            std::transform(order.data(), order.data() + order.size(),
+                    order.data(), [&](int index) { return adj_faces[index]; });
+
+            for (auto ei : adj_faces) {
+                const size_t fid = edge_index_to_face_index(ei);
+                patch_edge_adj[P[fid]].push_back(ei);
+            }
+        }
+    }
+
+    const int INVALID = std::numeric_limits<int>::max();
+    cells.resize(num_patches, 2);
+    cells.setConstant(INVALID);
+
+    auto peel_cell_bd = [&](
+            size_t seed_patch_id,
+            short seed_patch_side,
+            size_t cell_idx) {
+        typedef std::pair<size_t, short> PatchSide;
+        std::queue<PatchSide> Q;
+        Q.emplace(seed_patch_id, seed_patch_side);
+        cells(seed_patch_id, seed_patch_side) = cell_idx;
+        while (!Q.empty()) {
+            auto entry = Q.front();
+            Q.pop();
+            const size_t patch_id = entry.first;
+            const short side = entry.second;
+            const auto& adj_edges = patch_edge_adj[patch_id];
+            for (const auto& ei : adj_edges) {
+                const size_t uei = EMAP[ei];
+                const auto& order = orders[uei];
+                const auto& orientation = orientations[uei];
+                const size_t edge_valance = order.size();
+                size_t curr_i = 0;
+                for (curr_i=0; curr_i < edge_valance; curr_i++) {
+                    if (order[curr_i] == ei) break;
+                }
+                assert(curr_i < edge_valance);
+
+                const bool cons = orientation[curr_i];
+                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;
+                }
+                const size_t next_ei = order[next];
+                const bool next_cons = orientation[next];
+                const size_t next_patch_id = P[next_ei % num_faces];
+                const short next_patch_side = (next_cons != cons) ?
+                    side:abs(side-1);
+                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(cells(next_patch_id, next_patch_side) == cell_idx);
+                }
+            }
+        }
+    };
+
+    size_t count=0;
+    for (size_t i=0; i<num_patches; i++) {
+        if (cells(i, 0) == INVALID) {
+            peel_cell_bd(i, 0, count);
+            count++;
+        }
+        if (cells(i, 1) == INVALID) {
+            peel_cell_bd(i, 1, count);
+            count++;
+        }
+    }
+    return count;
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedE,
+    typename DeriveduE,
+    typename uE2EType,
+    typename DerivedEMAP,
+    typename DerivedC >
+IGL_INLINE size_t igl::cgal::extract_cells(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const Eigen::PlainObjectBase<DerivedE>& E,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        Eigen::PlainObjectBase<DerivedC>& cells) {
+    const size_t num_faces = F.rows();
+    typedef typename DerivedF::Scalar Index;
+    const size_t num_patches = P.maxCoeff()+1;
+
+    DerivedC raw_cells;
+    const size_t num_raw_cells =
+        igl::cgal::extract_cells_single_component(
+                V, F, P, uE, uE2E, EMAP, raw_cells);
+
+    std::vector<std::vector<std::vector<Index > > > TT,_1;
+    igl::triangle_triangle_adjacency(E, EMAP, uE2E, false, TT, _1);
+
+    Eigen::VectorXi C, counts;
+    igl::facet_components(TT, C, counts);
+
+    const size_t num_components = counts.size();
+    std::vector<std::vector<size_t> > components(num_components);
+    for (size_t i=0; i<num_faces; i++) {
+        components[C[i]].push_back(i);
+    }
+    std::vector<Eigen::VectorXi> 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());
+    }
+
+    Eigen::VectorXi outer_facets(num_components);
+    Eigen::VectorXi outer_facet_orientation(num_components);
+    Eigen::VectorXi outer_cells(num_components);
+    for (size_t i=0; i<num_components; i++) {
+        bool flipped;
+        igl::cgal::outer_facet(V, F, Is[i], outer_facets[i], flipped);
+        outer_facet_orientation[i] = flipped?1:0;
+        outer_cells[i] = raw_cells(P[outer_facets[i]], outer_facet_orientation[i]);
+    }
+
+    auto get_triangle_center = [&](const size_t fid) {
+        return ((V.row(F(fid, 0)) + V.row(F(fid, 1)) + V.row(F(fid, 2)))
+                /3.0).eval();
+    };
+    std::vector<std::vector<size_t> > nested_cells(num_raw_cells);
+    std::vector<std::vector<size_t> > ambient_cells(num_raw_cells);
+    std::vector<std::vector<size_t> > ambient_comps(num_components);
+    if (num_components > 1) {
+        for (size_t i=0; i<num_components; i++) {
+            DerivedV queries(num_components-1, 3);
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                size_t index = j<i ? j:j-1;
+                queries.row(index) = get_triangle_center(outer_facets[j]);
+            }
+
+            const auto& I = Is[i];
+            Eigen::VectorXi closest_facets, closest_facet_orientations;
+            igl::cgal::closest_facet(V, F, I, queries,
+                    closest_facets, closest_facet_orientations);
+
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                size_t index = j<i ? j:j-1;
+                const size_t closest_patch = P[closest_facets[index]];
+                const size_t closest_patch_side = closest_facet_orientations[index]
+                    ? 0:1;
+                const size_t ambient_cell = raw_cells(closest_patch,
+                        closest_patch_side);
+                if (ambient_cell != outer_cells[i]) {
+                    nested_cells[ambient_cell].push_back(outer_cells[j]);
+                    ambient_cells[outer_cells[j]].push_back(ambient_cell);
+                    ambient_comps[j].push_back(i);
+                }
+            }
+        }
+    }
+
+    const size_t INVALID = std::numeric_limits<size_t>::max();
+    const size_t INFINITE_CELL = num_raw_cells;
+    std::vector<size_t> embedded_cells(num_raw_cells, INVALID);
+    for (size_t i=0; i<num_components; i++) {
+        const size_t outer_cell = outer_cells[i];
+        const auto& ambient_comps_i = ambient_comps[i];
+        const auto& ambient_cells_i = ambient_cells[outer_cell];
+        const size_t num_ambient_comps = ambient_comps_i.size();
+        assert(num_ambient_comps == ambient_cells_i.size());
+        if (num_ambient_comps > 0) {
+            size_t embedded_comp = INVALID;
+            size_t embedded_cell = INVALID;
+            for (size_t j=0; j<num_ambient_comps; j++) {
+                if (ambient_comps[ambient_comps_i[j]].size() ==
+                        num_ambient_comps-1) {
+                    embedded_comp = ambient_comps_i[j];
+                    embedded_cell = ambient_cells_i[j];
+                    break;
+                }
+            }
+            assert(embedded_comp != INVALID);
+            assert(embedded_cell != INVALID);
+            embedded_cells[outer_cell] = embedded_cell;
+        } else {
+            embedded_cells[outer_cell] = INFINITE_CELL;
+        }
+    }
+    for (size_t i=0; i<num_patches; i++) {
+        if (embedded_cells[raw_cells(i,0)] != INVALID) {
+            raw_cells(i,0) = embedded_cells[raw_cells(i, 0)];
+        }
+        if (embedded_cells[raw_cells(i,1)] != INVALID) {
+            raw_cells(i,1) = embedded_cells[raw_cells(i, 1)];
+        }
+    }
+
+    size_t count = 0;
+    std::vector<size_t> mapped_indices(num_raw_cells+1, INVALID);
+    for (size_t i=0; i<num_patches; i++) {
+        const size_t old_positive_cell_id = raw_cells(i, 0);
+        const size_t old_negative_cell_id = raw_cells(i, 1);
+        size_t positive_cell_id, negative_cell_id;
+        if (mapped_indices[old_positive_cell_id] == INVALID) {
+            mapped_indices[old_positive_cell_id] = count;
+            positive_cell_id = count;
+            count++;
+        } else {
+            positive_cell_id = mapped_indices[old_positive_cell_id];
+        }
+        if (mapped_indices[old_negative_cell_id] == INVALID) {
+            mapped_indices[old_negative_cell_id] = count;
+            negative_cell_id = count;
+            count++;
+        } else {
+            negative_cell_id = mapped_indices[old_negative_cell_id];
+        }
+        raw_cells(i, 0) = positive_cell_id;
+        raw_cells(i, 1) = negative_cell_id;
+    }
+    cells = raw_cells;
+    return count;
+}
+
+template<
+typename DerivedV,
+typename DerivedF,
+typename DerivedC >
+IGL_INLINE size_t igl::cgal::extract_cells(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        Eigen::PlainObjectBase<DerivedC>& cells) {
+    const size_t num_faces = F.rows();
+    typedef typename DerivedF::Scalar Index;
+
+    Eigen::MatrixXi E, uE;
+    Eigen::VectorXi EMAP;
+    std::vector<std::vector<size_t> > uE2E;
+    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+    Eigen::VectorXi P;
+    const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+
+    DerivedC per_patch_cells;
+    const size_t num_cells =
+        igl::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+
+    cells.resize(num_faces, 2);
+    for (size_t i=0; i<num_faces; i++) {
+        cells.row(i) = per_patch_cells.row(P[i]);
+    }
+    return num_cells;
+}

include/igl/cgal/extract_cells.h

@@ -0,0 +1,111 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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_CGAL_EXTRACT_CELLS
+#define IGL_CGAL_EXTRACT_CELLS
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl {
+    namespace cgal {
+        // Extract connected 3D space partitioned by mesh (V, F).
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //
+        // Output:
+        //   cells  #F by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of face i, and cells(i,1)
+        //          represents cell index of the negqtive side.
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedC >
+        IGL_INLINE size_t extract_cells(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                Eigen::PlainObjectBase<DerivedC>& cells);
+
+        // Extract connected 3D space partitioned by mesh (V, F).
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //   P  #F list of patch indices.
+        //   E  #E by 2 array of vertex indices, one edge per row.
+        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+        //  EMAP  #F*3 list of indices into uE.
+        //
+        // Output:
+        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of patch i, and cells(i,1)
+        //          represents cell index of the negqtive side.
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedP,
+            typename DerivedE,
+            typename DeriveduE,
+            typename uE2EType,
+            typename DerivedEMAP,
+            typename DerivedC >
+        IGL_INLINE size_t extract_cells(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                const Eigen::PlainObjectBase<DerivedP>& P,
+                const Eigen::PlainObjectBase<DerivedE>& E,
+                const Eigen::PlainObjectBase<DeriveduE>& uE,
+                const std::vector<std::vector<uE2EType> >& uE2E,
+                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+                Eigen::PlainObjectBase<DerivedC>& cells);
+
+        // Extract connected 3D space partitioned by mesh (V, F).  Each
+        // connected component of the mesh partitions the ambient space
+        // separately.
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //   P  #F list of patch indices.
+        //   E  #E by 2 array of vertex indices, one edge per row.
+        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+        //  EMAP  #F*3 list of indices into uE.
+        //
+        // Output:
+        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of patch i, and cells(i,1)
+        //          represents cell index of the negqtive side.
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedP,
+            typename DeriveduE,
+            typename uE2EType,
+            typename DerivedEMAP,
+            typename DerivedC >
+        IGL_INLINE size_t extract_cells_single_component(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                const Eigen::PlainObjectBase<DerivedP>& P,
+                const Eigen::PlainObjectBase<DeriveduE>& uE,
+                const std::vector<std::vector<uE2EType> >& uE2E,
+                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+                Eigen::PlainObjectBase<DerivedC>& cells);
+
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "extract_cells.cpp"
+#endif
+#endif

include/igl/cgal/order_facets_around_edge.cpp

@@ -1,6 +1,8 @@
 #include "order_facets_around_edge.h"
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
+#include <stdexcept>
+
 namespace igl
 {
     namespace cgal
@@ -38,7 +40,7 @@ void igl::cgal::order_facets_around_edge(
     size_t s,
     size_t d, 
     const std::vector<int>& adj_faces,
-    Eigen::PlainObjectBase<DerivedI>& order)
+    Eigen::PlainObjectBase<DerivedI>& order, bool debug)
 {
     using namespace igl::cgal::order_facets_around_edges_helper;
 
@@ -131,7 +133,10 @@ void igl::cgal::order_facets_around_edge(
     const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
     const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
     const Plane_3 separator(p_s, p_d, p_o);
-    assert(!separator.is_degenerate());
+    if (separator.is_degenerate()) {
+        throw std::runtime_error(
+                "Cannot order facets around edge due to degenerated facets");
+    }
 
     std::vector<Point_3> opposite_vertices;
     for (size_t i=0; i<num_adj_faces; i++)
@@ -180,20 +185,43 @@ void igl::cgal::order_facets_around_edge(
                             break;
                         case CGAL::COLLINEAR:
                         default:
-                            assert(false);
+                            throw std::runtime_error(
+                                    "Degenerated facet detected.");
                             break;
                     }
                 }
                 break;
             default:
                 // Should not be here.
-                assert(false);
+                throw std::runtime_error("Unknown CGAL state detected.");
+        }
+    }
+    if (debug) {
+        std::cout << "tie positive: " << std::endl;
+        for (auto& f : tie_positive_oriented) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "positive side: " << std::endl;
+        for (auto& f : positive_side) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "tie negative: " << std::endl;
+        for (auto& f : tie_negative_oriented) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "negative side: " << std::endl;
+        for (auto& f : negative_side) {
+            std::cout << get_face_index(f) << " ";
         }
+        std::cout << std::endl;
     }
 
     Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
-    order_facets_around_edge(V, F, s, d, positive_side, positive_order);
-    order_facets_around_edge(V, F, s, d, negative_side, negative_order);
+    order_facets_around_edge(V, F, s, d, positive_side, positive_order, debug);
+    order_facets_around_edge(V, F, s, d, negative_side, negative_order, debug);
     std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
     std::vector<size_t> tie_negative_order = index_sort(tie_negative_oriented);
 
@@ -281,6 +309,7 @@ void igl::cgal::order_facets_around_edge(
 
     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;
@@ -303,7 +332,8 @@ void igl::cgal::order_facets_around_edge(
         K::Point_3 pd(V(d,0), V(d,1), V(d,2));
         K::Point_3 pp(pivot_point(0,0), pivot_point(0,1), pivot_point(0,2));
         if (CGAL::collinear(ps, pd, pp)) {
-            throw "Pivot point is collinear with the outer edge!";
+            throw std::runtime_error(
+                    "Pivot point is collinear with the outer edge!");
         }
     }
 
@@ -313,14 +343,16 @@ void igl::cgal::order_facets_around_edge(
     // Because face indices are used for tie breaking, the original face indices
     // in the new faces array must be ascending.
     auto comp = [&](int i, int j) {
-        return signed_index_to_index(i) < signed_index_to_index(j);
+        return signed_index_to_index(adj_faces[i]) <
+            signed_index_to_index(adj_faces[j]);
     };
-    std::vector<int> ordered_adj_faces(adj_faces);
-    std::sort(ordered_adj_faces.begin(), ordered_adj_faces.end(), comp);
+    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(ordered_adj_faces[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;
@@ -330,7 +362,7 @@ void igl::cgal::order_facets_around_edge(
     DerivedF faces(num_faces, 3);
     for (size_t i=0; i<N; i++)
     {
-        if (ordered_adj_faces[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  ;
@@ -377,15 +409,16 @@ void igl::cgal::order_facets_around_edge(
 
     for (size_t i=0; i<N; i++)
     {
-        order[i] = order_with_pivot[(pivot_index+i+1)%num_faces];
+        order[i] = adj_order[order_with_pivot[(pivot_index+i+1)%num_faces]];
     }
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::cgal::order_facets_around_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::order_facets_around_edge<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void
+igl::cgal::order_facets_around_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool);
+template void igl::cgal::order_facets_around_edge<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool);
 template void igl::cgal::order_facets_around_edge<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&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::cgal::order_facets_around_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::cgal::order_facets_around_edge<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/cgal/order_facets_around_edge.h

@@ -43,7 +43,8 @@ namespace igl {
                 const Eigen::PlainObjectBase<DerivedF>& F,
                 size_t s, size_t d, 
                 const std::vector<int>& adj_faces,
-                Eigen::PlainObjectBase<DerivedI>& order);
+                Eigen::PlainObjectBase<DerivedI>& order,
+                bool debug=false);
 
         // This funciton is a wrapper around the one above.  Since the ordering
         // is circular, the pivot point is used to define a starting point.  So

include/igl/cgal/peel_winding_number_layers.cpp

@@ -0,0 +1,33 @@
+#include "peel_winding_number_layers.h"
+
+#include <cassert>
+
+#include "propagate_winding_numbers.h"
+
+template<
+typename DerivedV,
+typename DerivedF,
+typename DerivedW >
+IGL_INLINE size_t igl::cgal::peel_winding_number_layers(
+        const Eigen::PlainObjectBase<DerivedV > & V,
+        const Eigen::PlainObjectBase<DerivedF > & F,
+        Eigen::PlainObjectBase<DerivedW>& W) {
+    const size_t num_faces = F.rows();
+    Eigen::VectorXi labels(num_faces);
+    labels.setZero();
+
+    Eigen::MatrixXi winding_numbers;
+    igl::cgal::propagate_winding_numbers(V, F, labels, winding_numbers);
+    assert(winding_numbers.rows() == num_faces);
+    assert(winding_numbers.cols() == 2);
+
+    int min_w = winding_numbers.minCoeff();
+    int max_w = winding_numbers.maxCoeff();
+    assert(max_w > min_w);
+
+    W.resize(num_faces, 1);
+    for (size_t i=0; i<num_faces; i++) {
+        W(i, 0) = winding_numbers(i, 1);
+    }
+    return max_w - min_w;
+}

include/igl/cgal/peel_winding_number_layers.h

@@ -0,0 +1,22 @@
+#ifndef IGL_CGAL_PEEL_WINDING_NUMBER_LAYERS_H
+#define IGL_CGAL_PEEL_WINDING_NUMBER_LAYERS_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+    namespace cgal {
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedW >
+        IGL_INLINE size_t peel_winding_number_layers(
+                const Eigen::PlainObjectBase<DerivedV > & V,
+                const Eigen::PlainObjectBase<DerivedF > & F,
+                Eigen::PlainObjectBase<DerivedW>& W);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "peel_winding_number_layers.cpp"
+#endif
+#endif

include/igl/cgal/propagate_winding_numbers.cpp

@@ -0,0 +1,257 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "propagate_winding_numbers.h"
+#include "../extract_manifold_patches.h"
+#include "../extract_non_manifold_edge_curves.h"
+#include "../facet_components.h"
+#include "../unique_edge_map.h"
+#include "../writeOBJ.h"
+#include "../writePLY.h"
+#include "order_facets_around_edge.h"
+#include "outer_facet.h"
+#include "closest_facet.h"
+#include "assign_scalar.h"
+#include "extract_cells.h"
+
+#include <stdexcept>
+#include <limits>
+#include <vector>
+#include <tuple>
+#include <queue>
+
+namespace propagate_winding_numbers_helper {
+    template<
+        typename DerivedF,
+        typename DeriveduE,
+        typename uE2EType >
+    bool is_orientable(
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DeriveduE>& uE,
+            const std::vector<std::vector<uE2EType> >& uE2E) {
+        const size_t num_faces = F.rows();
+        const size_t num_edges = uE.rows();
+        auto edge_index_to_face_index = [&](size_t ei) {
+            return ei % num_faces;
+        };
+        auto is_consistent = [&](size_t fid, size_t s, size_t d) {
+            if (F(fid, 0) == s && F(fid, 1) == d) return true;
+            if (F(fid, 1) == s && F(fid, 2) == d) return true;
+            if (F(fid, 2) == s && F(fid, 0) == d) return true;
+
+            if (F(fid, 0) == d && F(fid, 1) == s) return false;
+            if (F(fid, 1) == d && F(fid, 2) == s) return false;
+            if (F(fid, 2) == d && F(fid, 0) == s) return false;
+            throw "Invalid face!!";
+        };
+        for (size_t i=0; i<num_edges; i++) {
+            const size_t s = uE(i,0);
+            const size_t d = uE(i,1);
+            int count=0;
+            for (const auto& ei : uE2E[i]) {
+                const size_t fid = edge_index_to_face_index(ei);
+                if (is_consistent(fid, s, d)) count++;
+                else count--;
+            }
+            if (count != 0) {
+                return false;
+            }
+        }
+        return true;
+    }
+}
+
+template<
+typename DerivedV,
+typename DerivedF,
+typename DerivedL,
+typename DerivedW>
+IGL_INLINE void igl::cgal::propagate_winding_numbers(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedL>& labels,
+        Eigen::PlainObjectBase<DerivedW>& W) {
+    const size_t num_faces = F.rows();
+    typedef typename DerivedF::Scalar Index;
+
+    Eigen::MatrixXi E, uE;
+    Eigen::VectorXi EMAP;
+    std::vector<std::vector<size_t> > uE2E;
+    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+    assert(propagate_winding_numbers_helper::is_orientable(F, uE, uE2E));
+
+    Eigen::VectorXi P;
+    const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+
+    DerivedW per_patch_cells;
+    const size_t num_cells =
+        igl::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+
+    typedef std::tuple<size_t, bool, size_t> CellConnection;
+    std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+    for (size_t i=0; i<num_patches; i++) {
+        const int positive_cell = per_patch_cells(i,0);
+        const int negative_cell = per_patch_cells(i,1);
+        cell_adjacency[positive_cell].emplace(negative_cell, false, i);
+        cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+    }
+
+    auto save_cell = [&](const std::string& filename, size_t cell_id) {
+        std::vector<size_t> faces;
+        for (size_t i=0; i<num_patches; i++) {
+            if ((per_patch_cells.row(i).array() == cell_id).any()) {
+                for (size_t j=0; j<num_faces; j++) {
+                    if (P[j] == i) {
+                        faces.push_back(j);
+                    }
+                }
+            }
+        }
+        Eigen::MatrixXi cell_faces(faces.size(), 3);
+        for (size_t i=0; i<faces.size(); i++) {
+            cell_faces.row(i) = F.row(faces[i]);
+        }
+        Eigen::MatrixXd vertices(V.rows(), 3);
+        for (size_t i=0; i<V.rows(); i++) {
+            assign_scalar(V(i,0), vertices(i,0));
+            assign_scalar(V(i,1), vertices(i,1));
+            assign_scalar(V(i,2), vertices(i,2));
+        }
+        writePLY(filename, vertices, cell_faces);
+    };
+
+    {
+        // Check for odd cycle.
+        Eigen::VectorXi cell_labels(num_cells);
+        cell_labels.setZero();
+        Eigen::VectorXi parents(num_cells);
+        parents.setConstant(-1);
+        auto trace_parents = [&](size_t idx) {
+            std::list<size_t> path;
+            path.push_back(idx);
+            while (parents[path.back()] != path.back()) {
+                path.push_back(parents[path.back()]);
+            }
+            return path;
+        };
+        for (size_t i=0; i<num_cells; i++) {
+            if (cell_labels[i] == 0) {
+                cell_labels[i] = 1;
+                std::queue<size_t> Q;
+                Q.push(i);
+                parents[i] = i;
+                while (!Q.empty()) {
+                    size_t curr_idx = Q.front();
+                    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) {
+                            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) {
+                                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) {
+                                    std::cout << cell_id << " ";
+                                    std::stringstream filename;
+                                    filename << "cell_" << cell_id << ".ply";
+                                    save_cell(filename.str(), cell_id);
+                                }
+                                std::cout << std::endl;
+                            }
+                            assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    size_t outer_facet;
+    bool flipped;
+    Eigen::VectorXi I;
+    I.setLinSpaced(num_faces, 0, num_faces-1);
+    igl::cgal::outer_facet(V, F, I, outer_facet, flipped);
+
+    const size_t outer_patch = P[outer_facet];
+    const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
+
+    Eigen::VectorXi patch_labels(num_patches);
+    const int INVALID = std::numeric_limits<int>::max();
+    patch_labels.setConstant(INVALID);
+    for (size_t i=0; i<num_faces; i++) {
+        if (patch_labels[P[i]] == INVALID) {
+            patch_labels[P[i]] = labels[i];
+        } else {
+            assert(patch_labels[P[i]] == labels[i]);
+        }
+    }
+    assert((patch_labels.array() != INVALID).all());
+    const size_t num_labels = patch_labels.maxCoeff()+1;
+
+    Eigen::MatrixXi per_cell_W(num_cells, num_labels);
+    per_cell_W.setConstant(INVALID);
+    per_cell_W.row(infinity_cell).setZero();
+    std::queue<size_t> Q;
+    Q.push(infinity_cell);
+    while (!Q.empty()) {
+        size_t curr_cell = Q.front();
+        Q.pop();
+        for (const auto& neighbor : cell_adjacency[curr_cell]) {
+            size_t neighbor_cell, patch_idx;
+            bool direction;
+            std::tie(neighbor_cell, direction, patch_idx) = neighbor;
+            if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
+                per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
+                for (size_t i=0; i<num_labels; i++) {
+                    int inc = (patch_labels[patch_idx] == i) ?
+                        (direction ? -1:1) :0;
+                    per_cell_W(neighbor_cell, i) =
+                        per_cell_W(curr_cell, i) + inc;
+                }
+                Q.push(neighbor_cell);
+            } else {
+                for (size_t i=0; i<num_labels; i++) {
+                    if (i == patch_labels[patch_idx]) {
+                        int inc = direction ? -1:1;
+                        assert(per_cell_W(neighbor_cell, i) ==
+                                per_cell_W(curr_cell, i) + inc);
+                    } else {
+                        assert(per_cell_W(neighbor_cell, i) ==
+                                per_cell_W(curr_cell, i));
+                    }
+                }
+            }
+        }
+    }
+
+    W.resize(num_faces, num_labels*2);
+    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);
+        for (size_t j=0; j<num_labels; j++) {
+            W(i,j*2  ) = per_cell_W(positive_cell, j);
+            W(i,j*2+1) = per_cell_W(negative_cell, j);
+        }
+    }
+
+    //for (size_t i=0; i<num_cells; i++) {
+    //    std::stringstream filename;
+    //    filename << "cell_" << i << "_w_" << per_cell_W(i, 1) << ".ply";
+    //    save_cell(filename.str(), i);
+    //}
+}
+

include/igl/cgal/propagate_winding_numbers.h

@@ -0,0 +1,57 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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_CGAL_PROPAGATE_WINDING_NUMBERS_H
+#define IGL_CGAL_PROPAGATE_WINDING_NUMBERS_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+// The following methods compute the winding number on each side of each facet
+// or patch of a 3D mesh.  The input mesh is valid if it splits the ambient
+// space, R^3, into subspaces with constant integer winding numbers.  That is
+// the input mesh should be closed and for each directed edge the number of
+// clockwise facing facets should equal the number of counterclockwise facing
+// facets.
+
+namespace igl {
+    namespace cgal {
+
+        // 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
+        // closed, consistently oriented and induces integer winding numbers.
+        //
+        // Inputs:
+        //   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``.
+        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);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "propagate_winding_numbers.cpp"
+#endif
+#endif

include/igl/cgal/remesh_intersections.cpp

@@ -1,328 +1,338 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "remesh_intersections.h"
-#include "SelfIntersectMesh.h"
-#include "assign_scalar.h"
-#include "projected_delaunay.h"
+
+#include <vector>
+#include <map>
+#include <unordered_map>
 #include <iostream>
-#include <cassert>
 
 template <
-  typename DerivedV,
-  typename DerivedF,
-  typename Kernel,
-  typename DerivedVV,
-  typename DerivedFF,
-  typename DerivedJ,
-  typename DerivedIM>
+typename DerivedV,
+typename DerivedF,
+typename Kernel,
+typename DerivedVV,
+typename DerivedFF,
+typename DerivedJ,
+typename DerivedIM>
 IGL_INLINE void igl::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::pair<typename DerivedF::Index,
-      std::vector<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)
-{
-  using namespace std;
-  using namespace Eigen;
-  typedef typename DerivedF::Index          Index;
-  typedef CGAL::Point_3<Kernel>    Point_3;
-  //typedef CGAL::Segment_3<Kernel>  Segment_3; 
-  //typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-  typedef CGAL::Plane_3<Kernel>    Plane_3;
-  //typedef CGAL::Point_2<Kernel>    Point_2;
-  //typedef CGAL::Segment_2<Kernel>  Segment_2; 
-  //typedef CGAL::Triangle_2<Kernel> Triangle_2; 
-  typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
-  typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
-  typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
-  typedef CGAL::Exact_intersections_tag Itag;
-  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-    CDT_2;
-  typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
-  typedef std::pair<Index,Index> EMK;
-  typedef std::vector<Index> EMV;
-  //typedef std::map<EMK,EMV> EdgeMap;
-  typedef std::pair<Index,Index> EMK;
-  //typedef std::vector<CGAL::Object> ObjectList;
-  typedef std::vector<Index> IndexList;
+        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::pair<typename DerivedF::Index,
+        std::vector<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) {
+
+    typedef CGAL::Point_3<Kernel>    Point_3;
+    typedef CGAL::Segment_3<Kernel>  Segment_3; 
+    typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+    typedef CGAL::Plane_3<Kernel>    Plane_3;
+    typedef CGAL::Point_2<Kernel>    Point_2;
+    typedef CGAL::Segment_2<Kernel>  Segment_2; 
+    typedef CGAL::Triangle_2<Kernel> Triangle_2; 
+    typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+    typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
+    typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
+    typedef CGAL::Exact_intersections_tag Itag;
+    typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
+        CDT_2;
+    typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
+
+    typedef typename DerivedF::Index Index;
+    typedef std::pair<Index, Index> Edge;
+    struct EdgeHash {
+        size_t operator()(const Edge& e) const {
+            return (e.first * 805306457) ^ (e.second * 201326611);
+        }
+    };
+    typedef std::unordered_map<Edge, std::vector<Index>, EdgeHash > EdgeMap;
+
+    auto normalize_plane_coeff = [](const Plane_3& P) {
+        std::vector<typename Kernel::FT> coeffs = {
+            P.a(), P.b(), P.c(), P.d()
+        };
+        const auto max_itr = std::max_element(coeffs.begin(), coeffs.end());
+        const auto min_itr = std::min_element(coeffs.begin(), coeffs.end());
+        typename Kernel::FT max_coeff;
+        if (*max_itr < -1 * *min_itr) {
+            max_coeff = *min_itr;
+        } else {
+            max_coeff = *max_itr;
+        }
+        std::transform(coeffs.begin(), coeffs.end(), coeffs.begin(),
+                [&](const typename Kernel::FT& val)
+                {return val / max_coeff; } );
+        return coeffs;
+    };
+
+    auto plane_comp = [&](const Plane_3& p1, const Plane_3& p2) {
+        const auto p1_coeffs = normalize_plane_coeff(p1);
+        const auto p2_coeffs = normalize_plane_coeff(p2);
+        if (p1_coeffs[0] != p2_coeffs[0])
+            return p1_coeffs[0] < p2_coeffs[0];
+        if (p1_coeffs[1] != p2_coeffs[1])
+            return p1_coeffs[1] < p2_coeffs[1];
+        if (p1_coeffs[2] != p2_coeffs[2])
+            return p1_coeffs[2] < p2_coeffs[2];
+        if (p1_coeffs[3] != p2_coeffs[3])
+            return p1_coeffs[3] < p2_coeffs[3];
+        return false;
+    };
+    std::map<Plane_3, std::vector<Index>, decltype(plane_comp)>
+        unique_planes(plane_comp);
+
+    const size_t num_faces = F.rows();
+    const size_t num_base_vertices = V.rows();
+    assert(num_faces == T.size());
+    std::vector<bool> is_offending(num_faces, false);
+    for (const auto itr : offending) {
+        const auto& fid = itr.first;
+        is_offending[fid] = true;
 
-  int NF_count = 0;
-  // list of new faces, we'll fix F later
-  vector<
-    typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
-    > NF(offending.size());
-  // list of new vertices
-  typedef vector<Point_3> Point_3List;
-  Point_3List NV;
-  Index NV_count = 0;
-  vector<CDT_plus_2> cdt(offending.size());
-  vector<Plane_3> P(offending.size());
-  // Use map for *all* faces
-  map<typename CDT_plus_2::Vertex_handle,Index> v2i;
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  double t_proj_del = 0;
-#endif
-  // Unfortunately it looks like CGAL has trouble allocating memory when
-  // multiple openmp threads are running. Crashes durring CDT...
-  //// Loop over offending triangles
-  //const size_t noff = offending.size();
-//# pragma omp parallel for if (noff>1000)
-  for(const auto & okv : offending)
-  {
-    // index in F
-    const Index f = okv.first;
-    const Index o = okv.second.first;
-    {
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-      const double t_before = get_seconds();
-#endif
-      CDT_plus_2 cdt_o;
-      projected_delaunay(T[f],okv.second.second,cdt_o);
-      cdt[o] = cdt_o;
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-      t_proj_del += (get_seconds()-t_before);
-#endif
+        Plane_3 key = T[fid].supporting_plane();
+        assert(!key.is_degenerate());
+        const auto jtr = unique_planes.find(key);
+        if (jtr == unique_planes.end()) {
+            unique_planes.insert({key, {fid}});
+        } else {
+            jtr->second.push_back(fid);
+        }
     }
-    // Q: Is this also delaunay in 3D?
-    // A: No, because the projection is affine and delaunay is not affine
-    // invariant.
-    // Q: Then, can't we first get the 2D delaunay triangulation, then lift it
-    // to 3D and flip any offending edges?
-    // Plane of projection (also used by projected_delaunay)
-    P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2));
-    // Build index map
-    {
-      Index i=0;
-      for(
-        typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin();
-        vit != cdt[o].finite_vertices_end();
-        ++vit)
-      {
-        if(i<3)
-        {
-          //cout<<T[f].vertex(i)<<
-          //  (T[f].vertex(i) == P[o].to_3d(vit->point())?" == ":" != ")<<
-          //  P[o].to_3d(vit->point())<<endl;
-#ifndef NDEBUG
-          // I want to be sure that the original corners really show up as the
-          // original corners of the CDT. I.e. I don't trust CGAL to maintain
-          // the order
-          assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
-#endif
-          // For first three, use original index in F
-//#   pragma omp critical
-          v2i[vit] = F(f,i);
-        }else
-        {
-          const Point_3 vit_point_3 = P[o].to_3d(vit->point());
-          // First look up each edge's neighbors to see if exact point has
-          // already been added (This makes everything a bit quadratic)
-          bool found = false;
-          for(int e = 0; e<3 && !found;e++)
-          {
-            // Index of F's eth edge in V
-            Index i = F(f,(e+1)%3);
-            Index j = F(f,(e+2)%3);
-            // Be sure that i<j
-            if(i>j)
-            {
-              swap(i,j);
+
+    const size_t INVALID = std::numeric_limits<size_t>::max();
+    std::vector<std::vector<Index> > resolved_faces;
+    std::vector<Index> source_faces;
+    std::vector<Point_3> new_vertices;
+    EdgeMap edge_vertices;
+
+    /**
+     * Run constraint Delaunay triangulation on the plane.
+     */
+    auto run_delaunay_triangulation = [&](const Plane_3& P,
+            const std::vector<Index>& involved_faces,
+            std::vector<Point_3>& vertices,
+            std::vector<std::vector<Index> >& faces) {
+        CDT_plus_2 cdt;
+        for (const auto& fid : involved_faces) {
+            const auto itr = offending.find(fid);
+            const auto& triangle = T[fid];
+            cdt.insert_constraint(P.to_2d(triangle[0]), P.to_2d(triangle[1]));
+            cdt.insert_constraint(P.to_2d(triangle[1]), P.to_2d(triangle[2]));
+            cdt.insert_constraint(P.to_2d(triangle[2]), P.to_2d(triangle[0]));
+
+            if (itr == offending.end()) continue;
+            for (const auto& obj : itr->second.second) {
+                if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj)) {
+                    // Add segment constraint
+                    cdt.insert_constraint(
+                            P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+                }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj)) {
+                    // Add point
+                    cdt.insert(P.to_2d(*ipoint));
+                } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj)) {
+                    // Add 3 segment constraints
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(0)),P.to_2d(itri->vertex(1)));
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(1)),P.to_2d(itri->vertex(2)));
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(2)),P.to_2d(itri->vertex(0)));
+                } else if(const std::vector<Point_3 > *polyp = 
+                        CGAL::object_cast< std::vector<Point_3 > >(&obj)) {
+                    //cerr<<REDRUM("Poly...")<<endl;
+                    const std::vector<Point_3 > & poly = *polyp;
+                    const Index m = poly.size();
+                    assert(m>=2);
+                    for(Index p = 0;p<m;p++)
+                    {
+                        const Index np = (p+1)%m;
+                        cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
+                    }
+                }else {
+                    throw std::runtime_error("Unknown intersection object!");
+                }
             }
-            assert(edge2faces.count(EMK(i,j))==1);
-            const EMV & facesij = edge2faces.find(EMK(i,j))->second;
-            // loop over neighbors
-            for(
-              typename IndexList::const_iterator nit = facesij.begin();
-              nit != facesij.end() && !found;
-              nit++)
-            {
-              // don't consider self
-              if(*nit == f)
-              {
-                continue;
-              }
-              // index of neighbor in offending (to find its cdt)
-              assert(offending.count(*nit) == 1);
-              Index no = offending.find(*nit)->second.first;
-              // Loop over vertices of that neighbor's cdt (might not have been
-              // processed yet, but then it's OK because it'll just be empty)
-              for(
-                  typename CDT_plus_2::Finite_vertices_iterator uit = cdt[no].finite_vertices_begin();
-                  uit != cdt[no].finite_vertices_end() && !found;
-                  ++uit)
-              {
-                if(vit_point_3 == P[no].to_3d(uit->point()))
-                {
-                  assert(v2i.count(uit) == 1);
-//#   pragma omp critical
-                  v2i[vit] = v2i[uit];
-                  found = true;
+        }
+        const size_t num_vertices = cdt.number_of_vertices();
+        const size_t num_faces = cdt.number_of_faces();
+        std::map<typename CDT_plus_2::Vertex_handle,Index> v2i;
+        size_t count=0;
+        for (auto itr = cdt.finite_vertices_begin();
+                itr != cdt.finite_vertices_end(); itr++) {
+            vertices.push_back(P.to_3d(itr->point()));
+            v2i[itr] = count;
+            count++;
+        }
+        for (auto itr = cdt.finite_faces_begin();
+                itr != cdt.finite_faces_end(); itr++) {
+            faces.push_back( {
+                    v2i[itr->vertex(0)],
+                    v2i[itr->vertex(1)],
+                    v2i[itr->vertex(2)] });
+        }
+    };
+
+    /**
+     * 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 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];
+        }
+
+        // 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;
                 }
-              }
             }
-          }
-          if(!found)
-          {
-//#   pragma omp critical
-            {
-              v2i[vit] = V.rows()+NV_count;
-              NV.push_back(vit_point_3);
-              NV_count++;
+        }
+
+        // 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;
+    };
+
+    /**
+     * Determine the vertex indices for each corner of each output triangle.
+     */
+    auto post_triangulation_process = [&](
+            const std::vector<Point_3>& vertices,
+            const std::vector<std::vector<Index> >& faces,
+            const std::vector<Index>& involved_faces) {
+        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);
+                }
             }
-          }
         }
-        i++;
-      }
+    };
+
+    // Process un-touched faces.
+    for (size_t i=0; i<num_faces; i++) {
+        if (!is_offending[i]) {
+            resolved_faces.push_back(
+                    { F(i,0), F(i,1), F(i,2) } );
+            source_faces.push_back(i);
+        }
     }
-    {
-      Index i = 0;
-      // Resize to fit new number of triangles
-      NF[o].resize(cdt[o].number_of_faces(),3);
-//#   pragma omp atomic
-      NF_count+=NF[o].rows();
-      // Append new faces to NF
-      for(
-        typename CDT_plus_2::Finite_faces_iterator fit = cdt[o].finite_faces_begin();
-        fit != cdt[o].finite_faces_end();
-        ++fit)
-      {
-        NF[o](i,0) = v2i[fit->vertex(0)];
-        NF[o](i,1) = v2i[fit->vertex(1)];
-        NF[o](i,2) = v2i[fit->vertex(2)];
-        i++;
-      }
+
+    // Process self-intersecting faces.
+    for (const auto itr : unique_planes) {
+        Plane_3 P = itr.first;
+        const auto& involved_faces = itr.second;
+
+        std::vector<Point_3> vertices;
+        std::vector<std::vector<Index> > faces;
+        run_delaunay_triangulation(P, involved_faces, vertices, faces);
+        post_triangulation_process(vertices, faces, involved_faces);
     }
-  }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
-#endif
 
-  assert(NV_count == (Index)NV.size());
-  // Build output
-#ifndef NDEBUG
-  //{
-  //  Index off_count = 0;
-  //  for(Index f = 0;f<F.rows();f++)
-  //  {
-  //    off_count+= (offensive[f]?1:0);
-  //  }
-  //  assert(off_count==(Index)offending.size());
-  //}
-#endif
-  // Append faces
-  FF.resize(F.rows()-offending.size()+NF_count,3);
-  J.resize(FF.rows());
-  // First append non-offending original faces
-  // There's an Eigen way to do this in one line but I forget
-  Index off = 0;
-  for(Index f = 0;f<F.rows();f++)
-  {
-    if(!offending.count(f))
-    {
-      FF.row(off) = F.row(f);
-      J(off) = f;
-      off++;
+    // Output resolved mesh.
+    const size_t num_out_vertices = new_vertices.size() + num_base_vertices;
+    VV.resize(num_out_vertices, 3);
+    for (size_t i=0; i<num_base_vertices; i++) {
+        assign_scalar(V(i,0), VV(i,0));
+        assign_scalar(V(i,1), VV(i,1));
+        assign_scalar(V(i,2), VV(i,2));
     }
-  }
-  assert(off == (Index)(F.rows()-offending.size()));
-  // Now append replacement faces for offending faces
-  for(const auto & okv : offending)
-  {
-    // index in F
-    const Index f = okv.first;
-    const Index o = okv.second.first;
-    FF.block(off,0,NF[o].rows(),3) = NF[o];
-    J.block(off,0,NF[o].rows(),1).setConstant(f);
-    off += NF[o].rows();
-  }
-  // Append vertices
-  VV.resize(V.rows()+NV_count,3);
-  VV.block(0,0,V.rows(),3) = V.template cast<typename DerivedVV::Scalar>();
-  {
-    Index i = 0;
-    for(
-      typename Point_3List::const_iterator nvit = NV.begin();
-      nvit != NV.end();
-      nvit++)
-    {
-      for(Index d = 0;d<3;d++)
-      {
-        const Point_3 & p = *nvit;
-        // Don't convert via double if output type is same as Kernel
-        assign_scalar(p[d], VV(V.rows()+i,d));
-      }
-      i++;
+    for (size_t i=num_base_vertices; i<num_out_vertices; i++) {
+        assign_scalar(new_vertices[i-num_base_vertices][0], VV(i,0));
+        assign_scalar(new_vertices[i-num_base_vertices][1], VV(i,1));
+        assign_scalar(new_vertices[i-num_base_vertices][2], VV(i,2));
     }
-  }
-  IM.resize(VV.rows(),1);
-  map<Point_3,Index> vv2i;
-  // Safe to check for duplicates using double for original vertices: if
-  // incoming reps are different then the points are unique.
-  for(Index v = 0;v<V.rows();v++)
-  {
-    typename Kernel::FT p0,p1,p2;
-    assign_scalar(V(v,0),p0);
-    assign_scalar(V(v,1),p1);
-    assign_scalar(V(v,2),p2);
-    const Point_3 p(p0,p1,p2);
-    if(vv2i.count(p)==0)
-    {
-      vv2i[p] = v;
+
+    const size_t num_out_faces = resolved_faces.size();
+    FF.resize(num_out_faces, 3);
+    for (size_t i=0; i<num_out_faces; i++) {
+        FF(i,0) = resolved_faces[i][0];
+        FF(i,1) = resolved_faces[i][1];
+        FF(i,2) = resolved_faces[i][2];
     }
-    assert(vv2i.count(p) == 1);
-    IM(v) = vv2i[p];
-  }
-  // Must check for duplicates of new vertices using exact.
-  {
-    Index v = V.rows();
-    for(
-      typename Point_3List::const_iterator nvit = NV.begin();
-      nvit != NV.end();
-      nvit++)
-    {
-      const Point_3 & p = *nvit;
-      if(vv2i.count(p)==0)
-      {
-        vv2i[p] = v;
-      }
-      assert(vv2i.count(p) == 1);
-      IM(v) = vv2i[p];
-      v++;
+
+    J.resize(num_out_faces);
+    std::copy(source_faces.begin(), source_faces.end(), J.data());
+
+    // Extract unique vertex indices.
+    IM.resize(VV.rows(),1);
+    std::map<Point_3,Index> vv2i;
+    // Safe to check for duplicates using double for original vertices: if
+    // incoming reps are different then the points are unique.
+    for(Index v = 0;v<VV.rows();v++) {
+        typename Kernel::FT p0,p1,p2;
+        assign_scalar(VV(v,0),p0);
+        assign_scalar(VV(v,1),p1);
+        assign_scalar(VV(v,2),p2);
+        const Point_3 p(p0,p1,p2);
+        if(vv2i.count(p)==0) {
+            vv2i[p] = v;
+        }
+        assert(vv2i.count(p) == 1);
+        IM(v) = vv2i[p];
     }
-  }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"Output + dupes: "<<tictoc()<<endl;
-#endif
 }
 
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> const, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> const, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-#endif

include/igl/cgal/remesh_intersections.h

@@ -1,24 +1,18 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@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_CGAL_REMESH_INTERSECTIONS_H
 #define IGL_CGAL_REMESH_INTERSECTIONS_H
-#include <igl/igl_inline.h>
 
+#include <igl/igl_inline.h>
 #include <Eigen/Dense>
 #include "CGAL_includes.hpp"
 
-#ifdef MEX
-#  include <mex.h>
-#  include <cassert>
-#  undef assert
-#  define assert( isOK ) ( (isOK) ? (void)0 : (void) mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
-#endif
-  
 namespace igl
 {
   namespace cgal
@@ -72,6 +66,5 @@ namespace igl
 #ifndef IGL_STATIC_LIBRARY
 #  include "remesh_intersections.cpp"
 #endif
-  
-#endif
 
+#endif

include/igl/extract_manifold_patches.cpp

@@ -0,0 +1,68 @@
+#include "extract_manifold_patches.h"
+#include <cassert>
+#include <limits>
+#include <queue>
+
+template<
+typename DerivedF,
+typename DerivedEMAP,
+typename uE2EType,
+typename DerivedP>
+IGL_INLINE size_t igl::extract_manifold_patches(
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        Eigen::PlainObjectBase<DerivedP>& P) {
+    assert(F.cols() == 3);
+    const size_t num_faces = F.rows();
+
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto face_and_corner_index_to_edge_index = [&](size_t fi, size_t ci) {
+        return ci*num_faces + fi;
+    };
+    auto is_manifold_edge = [&](size_t fi, size_t ci) {
+        const size_t ei = face_and_corner_index_to_edge_index(fi, ci);
+        return uE2E[EMAP(ei, 0)].size() == 2;
+    };
+    auto get_adj_face_index = [&](size_t fi, size_t ci) -> size_t {
+        const size_t ei = face_and_corner_index_to_edge_index(fi, ci);
+        const auto& adj_faces = uE2E[EMAP(ei, 0)];
+        assert(adj_faces.size() == 2);
+        if (adj_faces[0] == ei) {
+            return edge_index_to_face_index(adj_faces[1]);
+        } else {
+            assert(adj_faces[1] == ei);
+            return edge_index_to_face_index(adj_faces[0]);
+        }
+    };
+
+    typedef typename DerivedP::Scalar Scalar;
+    const Scalar INVALID = std::numeric_limits<Scalar>::max();
+    P.resize(num_faces,1);
+    P.setConstant(INVALID);
+    size_t num_patches = 0;
+    for (size_t i=0; i<num_faces; i++) {
+        if (P(i,0) != INVALID) continue;
+
+        std::queue<size_t> Q;
+        Q.push(i);
+        P(i,0) = num_patches;
+        while (!Q.empty()) {
+            const size_t fid = Q.front();
+            Q.pop();
+            for (size_t j=0; j<3; j++) {
+                if (is_manifold_edge(fid, j)) {
+                    const size_t adj_fid = get_adj_face_index(fid, j);
+                    if (P(adj_fid,0) == INVALID) {
+                        Q.push(adj_fid);
+                        P(adj_fid,0) = num_patches;
+                    }
+                }
+            }
+        }
+        num_patches++;
+    }
+    assert((P.array() != INVALID).all());
+
+    return num_patches;
+}

include/igl/extract_manifold_patches.h

@@ -0,0 +1,42 @@
+#ifndef IGL_EXTRACT_MANIFOLD_PATCHES
+#define IGL_EXTRACT_MANIFOLD_PATCHES
+
+#include "igl_inline.h"
+#include <Eigen/Dense>
+#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.
+    //
+    // Inputs:
+    //   F  #F by 3 list representing triangles.
+    //   EMAP  #F*3 list of indices of unique undirected edges.
+    //   uE2E  #uE list of lists of indices into E of coexisting edges.
+    //
+    // Output:
+    //   P  #F list of patch incides.
+    //
+    // Returns:
+    //   number of manifold patches.
+    template <
+        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);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "extract_manifold_patches.cpp"
+#endif
+
+#endif

include/igl/extract_non_manifold_edge_curves.cpp

@@ -0,0 +1,116 @@
+#include "extract_non_manifold_edge_curves.h"
+#include <algorithm>
+#include <cassert>
+#include <list>
+#include <vector>
+#include <unordered_map>
+
+template<
+typename DerivedF,
+typename DerivedEMAP,
+typename uE2EType >
+IGL_INLINE void igl::extract_non_manifold_edge_curves(
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        std::vector<std::vector<size_t> >& curves) {
+    const size_t num_faces = F.rows();
+    assert(F.cols() == 3);
+    typedef std::pair<size_t, size_t> Edge;
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
+    auto get_edge_end_points = [&](size_t ei, size_t& s, size_t& d) {
+        const size_t fi = edge_index_to_face_index(ei);
+        const size_t ci = edge_index_to_corner_index(ei);
+        s = F(fi, (ci+1)%3);
+        d = F(fi, (ci+2)%3);
+    };
+
+    curves.clear();
+    const size_t num_unique_edges = uE2E.size();
+    std::unordered_multimap<size_t, size_t> vertex_edge_adjacency;
+    std::vector<size_t> non_manifold_edges;
+    for (size_t i=0; i<num_unique_edges; i++) {
+        const auto& adj_edges = uE2E[i];
+        if (adj_edges.size() == 2) continue;
+
+        const size_t ei = adj_edges[0];
+        size_t s,d;
+        get_edge_end_points(ei, s, d);
+        vertex_edge_adjacency.insert({{s, i}, {d, i}});
+        non_manifold_edges.push_back(i);
+        assert(vertex_edge_adjacency.count(s) > 0);
+        assert(vertex_edge_adjacency.count(d) > 0);
+    }
+
+    auto expand_forward = [&](std::list<size_t>& edge_curve,
+            size_t& front_vertex, size_t& end_vertex) {
+        while(vertex_edge_adjacency.count(front_vertex) == 2 &&
+                front_vertex != end_vertex) {
+            auto adj_edges = vertex_edge_adjacency.equal_range(front_vertex);
+            for (auto itr = adj_edges.first; itr!=adj_edges.second; itr++) {
+                const size_t uei = itr->second;
+                assert(uE2E.at(uei).size() != 2);
+                const size_t ei = uE2E[uei][0];
+                if (uei == edge_curve.back()) continue;
+                size_t s,d;
+                get_edge_end_points(ei, s, d);
+                edge_curve.push_back(uei);
+                if (s == front_vertex) {
+                    front_vertex = d;
+                } else if (d == front_vertex) {
+                    front_vertex = s;
+                } else {
+                    throw "Invalid vertex/edge adjacency!";
+                }
+                break;
+            }
+        }
+    };
+
+    auto expand_backward = [&](std::list<size_t>& edge_curve,
+            size_t& front_vertex, size_t& end_vertex) {
+        while(vertex_edge_adjacency.count(front_vertex) == 2 &&
+                front_vertex != end_vertex) {
+            auto adj_edges = vertex_edge_adjacency.equal_range(front_vertex);
+            for (auto itr = adj_edges.first; itr!=adj_edges.second; itr++) {
+                const size_t uei = itr->second;
+                assert(uE2E.at(uei).size() != 2);
+                const size_t ei = uE2E[uei][0];
+                if (uei == edge_curve.front()) continue;
+                size_t s,d;
+                get_edge_end_points(ei, s, d);
+                edge_curve.push_front(uei);
+                if (s == front_vertex) {
+                    front_vertex = d;
+                } else if (d == front_vertex) {
+                    front_vertex = s;
+                } else {
+                    throw "Invalid vertex/edge adjcency!";
+                }
+                break;
+            }
+        }
+    };
+
+    std::vector<bool> visited(num_unique_edges, false);
+    for (const size_t i : non_manifold_edges) {
+        if (visited[i]) continue;
+        std::list<size_t> edge_curve;
+        edge_curve.push_back(i);
+
+        const auto& adj_edges = uE2E[i];
+        assert(adj_edges.size() != 2);
+        const size_t ei = adj_edges[0];
+        size_t s,d;
+        get_edge_end_points(ei, s, d);
+
+        expand_forward(edge_curve, d, s);
+        expand_backward(edge_curve, s, d);
+        curves.emplace_back(edge_curve.begin(), edge_curve.end());
+        for (auto itr = edge_curve.begin(); itr!=edge_curve.end(); itr++) {
+            visited[*itr] = true;
+        }
+    }
+
+}

include/igl/extract_non_manifold_edge_curves.h

@@ -0,0 +1,40 @@
+#ifndef IGL_CGAL_EXTRACT_NON_MANIFOLD_EDGE_CURVES
+#define IGL_CGAL_EXTRACT_NON_MANIFOLD_EDGE_CURVES
+
+#include "igl_inline.h"
+#include <Eigen/Dense>
+#include <vector>
+
+namespace igl {
+    // Extract non-manifold curves from a given mesh.
+    // A non-manifold curves are a set of connected non-manifold edges that
+    // does not touch other non-manifold edges except at the end points.
+    // They are also maximal in the sense that they cannot be expanded by
+    // including more edges.
+    //
+    // Assumes the input mesh have all self-intersection resolved.  See
+    // ``igl::cgal::remesh_self_intersection`` for more details.
+    //
+    // Inputs:
+    //   F  #F by 3 list representing triangles.
+    //   EMAP  #F*3 list of indices of unique undirected edges.
+    //   uE2E  #uE list of lists of indices into E of coexisting edges.
+    //
+    // Output:
+    //   curves  An array of arries of unique edge indices.
+    template<
+        typename DerivedF,
+        typename DerivedEMAP,
+        typename uE2EType>
+    IGL_INLINE void extract_non_manifold_edge_curves(
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+            const std::vector<std::vector<uE2EType> >& uE2E,
+            std::vector<std::vector<size_t> >& curves);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "extract_non_manifold_edge_curves.cpp"
+#endif
+
+#endif