Merge branch 'qnzhou-master'

Former-commit-id: af249d8efc1849f2030098ab3d10a4fff940a715

include/igl/barycenter.h

@@ -11,17 +11,14 @@
 #include <Eigen/Dense>
 namespace igl
 {
-  // BARYCENTER
-  //
-  // B = barycenter(V,F)
-  //
-  // Compute the barycenter of every simplex
+  // Computes the barycenter of every simplex
   //
   // Inputs:
-  //   V #V x dim matrix of vertex coordinates
-  //   F #F x simplex_size  matrix of indices of simplex corners
+  //   V  #V x dim matrix of vertex coordinates
+  //   F  #F x simplex_size  matrix of indices of simplex corners into V
   // Output:
-  //   BC a #F x dim matrix of 3d vertices
+  //   BC  #F x dim matrix of 3d vertices
+  //
   template <
     typename DerivedV,
     typename DerivedF,

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,378 @@
 // 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/remesh_self_intersections.h>
-#include <igl/remove_unreferenced.h>
-#include <igl/mod.h>
-#include <igl/unique_simplices.h>
+#include "BinaryWindingNumberOperations.h"
+#include "../cgal/assign_scalar.h"
+#include "../cgal/propagate_winding_numbers.h"
+#include "../cgal/remesh_self_intersections.h"
+#include "../remove_unreferenced.h"
+#include "../unique_simplices.h"
+#include "../slice.h"
+
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <iostream>
+#include <algorithm>
 
-//#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;
-  }
 
-  // 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>&);
-
-
+template void igl::boolean::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#undef IGL_STATIC_LIBRARY
+#include "../remove_unreferenced.cpp"
+template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#include "../slice.cpp"
+template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > >(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&, int, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
 #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,360 @@
+// 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);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::cgal::closest_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

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/complex_to_mesh.cpp

@@ -140,8 +140,10 @@ IGL_INLINE bool igl::cgal::complex_to_mesh(
   }
 
   // CGAL code somehow can end up with unreferenced vertices
-  VectorXi _1;
-  remove_unreferenced( MatrixXd(V), MatrixXi(F), V,F,_1);
+  {
+    VectorXi _1;
+    remove_unreferenced( MatrixXd(V), MatrixXi(F), V,F,_1);
+  }
 
   return success;
 }

include/igl/cgal/extract_cells.cpp

@@ -0,0 +1,345 @@
+// 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;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+template unsigned long igl::cgal::extract_cells<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/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/mesh_to_cgal_triangle_list.cpp

@@ -60,4 +60,6 @@ template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, -1
 template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick>(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> > >&);
 template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck>(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> > >&);
 template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick>(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> > >&);
+template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick>(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&, std::__1::vector<CGAL::Triangle_3<CGAL::Epick>, std::__1::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+template void igl::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck>(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&, std::__1::vector<CGAL::Triangle_3<CGAL::Epeck>, std::__1::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
 #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,17 +409,18 @@ 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<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<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<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> >&);
 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<double, -1, -1, 0, -1, -1> > 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, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, unsigned long, unsigned long, std::__1::vector<int, std::__1::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

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/outer_element.cpp

@@ -0,0 +1,283 @@
+// 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 "outer_element.h"
+#include <iostream>
+#include <vector>
+
+template <
+     typename DerivedV,
+     typename DerivedF,
+     typename DerivedI,
+     typename IndexType,
+     typename DerivedA
+     >
+IGL_INLINE void igl::cgal::outer_vertex(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & v_index,
+        Eigen::PlainObjectBase<DerivedA> & A)
+{
+    // Algorithm: 
+    //    Find an outer vertex (i.e. vertex reachable from infinity)
+    //    Return the vertex with the largest X value.
+    //    If there is a tie, pick the one with largest Y value.
+    //    If there is still a tie, pick the one with the largest Z value.
+    //    If there is still a tie, then there are duplicated vertices within the
+    //    mesh, which violates the precondition.
+    typedef typename DerivedF::Scalar Index;
+    const Index INVALID = std::numeric_limits<Index>::max();
+    const size_t num_selected_faces = I.rows();
+    std::vector<size_t> candidate_faces;
+    Index outer_vid = INVALID;
+    typename DerivedV::Scalar outer_val = 0;
+    for (size_t i=0; i<num_selected_faces; i++)
+    {
+        size_t f = I(i);
+        for (size_t j=0; j<3; j++)
+        {
+            Index v = F(f, j);
+            auto vx = V(v, 0);
+            if (outer_vid == INVALID || vx > outer_val)
+            {
+                outer_val = vx;
+                outer_vid = v;
+                candidate_faces = {f};
+            } else if (v == outer_vid)
+            {
+                candidate_faces.push_back(f);
+            } else if (vx == outer_val)
+            {
+                // Break tie.
+                auto vy = V(v,1);
+                auto vz = V(v, 2);
+                auto outer_y = V(outer_vid, 1);
+                auto outer_z = V(outer_vid, 2);
+                assert(!(vy == outer_y && vz == outer_z));
+                bool replace = (vy > outer_y) ||
+                    ((vy == outer_y) && (vz > outer_z));
+                if (replace)
+                {
+                    outer_val = vx;
+                    outer_vid = v;
+                    candidate_faces = {f};
+                }
+            }
+        }
+    }
+
+    assert(outer_vid != INVALID);
+    assert(candidate_faces.size() > 0);
+    v_index = outer_vid;
+    A.resize(candidate_faces.size());
+    std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI,
+    typename IndexType,
+    typename DerivedA
+    >
+IGL_INLINE void igl::cgal::outer_edge(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & v1,
+        IndexType & v2,
+        Eigen::PlainObjectBase<DerivedA> & A) {
+    // Algorithm:
+    //    Find an outer vertex first.
+    //    Find the incident edge with largest abs slope when projected onto XY plane.
+    //    If there is a tie, check the signed slope and use the positive one.
+    //    If there is still a tie, break it using the projected slope onto ZX plane.
+    //    If there is still a tie, again check the signed slope and use the positive one.
+    //    If there is still a tie, then there are multiple overlapping edges,
+    //    which violates the precondition.
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
+    typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
+    const Index INVALID = std::numeric_limits<Index>::max();
+
+    Index outer_vid;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
+    outer_vertex(V, F, I, outer_vid, candidate_faces);
+    const ScalarArray3& outer_v = V.row(outer_vid);
+    assert(candidate_faces.size() > 0);
+
+    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
+    {
+        if (f[0] == vid) return 0;
+        if (f[1] == vid) return 1;
+        if (f[2] == vid) return 2;
+        assert(false);
+        return -1;
+    };
+
+    auto unsigned_value = [](Scalar v) -> Scalar {
+        if (v < 0) return v * -1;
+        else return v;
+    };
+
+    Scalar outer_slope_YX = 0;
+    Scalar outer_slope_ZX = 0;
+    Index outer_opp_vid = INVALID;
+    bool infinite_slope_detected = false;
+    std::vector<Index> incident_faces;
+    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
+        if (opp_vid == outer_opp_vid)
+        {
+            incident_faces.push_back(fid);
+            return;
+        }
+
+        const ScalarArray3 opp_v = V.row(opp_vid);
+        if (!infinite_slope_detected && outer_v[0] != opp_v[0])
+        {
+            // Finite slope
+            const ScalarArray3 diff = opp_v - outer_v;
+            const Scalar slope_YX = diff[1] / diff[0];
+            const Scalar slope_ZX = diff[2] / diff[0];
+            const Scalar u_slope_YX = unsigned_value(slope_YX);
+            const Scalar u_slope_ZX = unsigned_value(slope_ZX);
+            bool update = false;
+            if (outer_opp_vid == INVALID) {
+                update = true;
+            } else {
+                const Scalar u_outer_slope_YX = unsigned_value(outer_slope_YX);
+                if (u_slope_YX > u_outer_slope_YX) {
+                    update = true;
+                } else if (u_slope_YX == u_outer_slope_YX &&
+                        slope_YX > outer_slope_YX) {
+                    update = true;
+                } else if (slope_YX == outer_slope_YX) {
+                    const Scalar u_outer_slope_ZX =
+                        unsigned_value(outer_slope_ZX);
+                    if (u_slope_ZX > u_outer_slope_ZX) {
+                        update = true;
+                    } else if (u_slope_ZX == u_outer_slope_ZX &&
+                            slope_ZX > outer_slope_ZX) {
+                        update = true;
+                    } else if (slope_ZX == u_outer_slope_ZX) {
+                        assert(false);
+                    }
+                }
+            }
+
+            if (update) {
+                outer_opp_vid = opp_vid;
+                outer_slope_YX = slope_YX;
+                outer_slope_ZX = slope_ZX;
+                incident_faces = {fid};
+            }
+        } else if (!infinite_slope_detected)
+        {
+            // Infinite slope
+            outer_opp_vid = opp_vid;
+            infinite_slope_detected = true;
+            incident_faces = {fid};
+        }
+    };
+
+    const size_t num_candidate_faces = candidate_faces.size();
+    for (size_t i=0; i<num_candidate_faces; i++)
+    {
+        const Index fid = candidate_faces(i);
+        const IndexArray3& f = F.row(fid);
+        size_t id = get_vertex_index(f, outer_vid);
+        Index next_vid = f((id+1)%3);
+        Index prev_vid = f((id+2)%3);
+        check_and_update_outer_edge(next_vid, fid);
+        check_and_update_outer_edge(prev_vid, fid);
+    }
+
+    v1 = outer_vid;
+    v2 = outer_opp_vid;
+    A.resize(incident_faces.size());
+    std::copy(incident_faces.begin(), incident_faces.end(), A.data());
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedI,
+    typename IndexType
+    >
+IGL_INLINE void igl::cgal::outer_facet(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedN> & N,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & f,
+        bool & flipped) {
+    // Algorithm:
+    //    Find an outer edge.
+    //    Find the incident facet with the largest absolute X normal component.
+    //    If there is a tie, keep the one with positive X component.
+    //    If there is still a tie, pick the face with the larger signed index
+    //    (flipped face has negative index).
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    const size_t INVALID = std::numeric_limits<size_t>::max();
+
+    Index v1,v2;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;
+    outer_edge(V, F, I, v1, v2, incident_faces);
+    assert(incident_faces.size() > 0);
+
+    auto generic_fabs = [&](const Scalar& val) -> const Scalar {
+        if (val >= 0) return val;
+        else return -val;
+    };
+
+    Scalar max_nx = 0;
+    size_t outer_fid = INVALID;
+    const size_t num_incident_faces = incident_faces.size();
+    for (size_t i=0; i<num_incident_faces; i++) 
+    {
+        const auto& fid = incident_faces(i);
+        const Scalar nx = N(fid, 0);
+        if (outer_fid == INVALID) {
+            max_nx = nx;
+            outer_fid = fid;
+        } else {
+            if (generic_fabs(nx) > generic_fabs(max_nx)) {
+                max_nx = nx;
+                outer_fid = fid;
+            } else if (nx == -max_nx && nx > 0) {
+                max_nx = nx;
+                outer_fid = fid;
+            } else if (nx == max_nx) {
+                if ((max_nx >= 0 && outer_fid < fid) ||
+                    (max_nx <  0 && outer_fid > fid)) {
+                    max_nx = nx;
+                    outer_fid = fid;
+                }
+            }
+        }
+    }
+
+    assert(outer_fid != INVALID);
+    f = outer_fid;
+    flipped = max_nx < 0;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+template void igl::cgal::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, unsigned long&, bool&);
+template void igl::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::cgal::outer_facet<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>, int>(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&, int&, bool&);
+template void igl::cgal::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> >&);
+template void igl::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/cgal/outer_element.h

@@ -0,0 +1,112 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingan 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_OUTER_ELEMENT_H
+#define IGL_CGAL_OUTER_ELEMENT_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  namespace cgal
+  {
+    // Find a vertex that is reachable from infinite without crossing any faces.
+    // Such vertex is called "outer vertex."
+    //
+    // Precondition: The input mesh must have all self-intersection resolved and
+    // no duplicated vertices.  See cgal::remesh_self_intersections.h for how to
+    // obtain such input.
+    //
+    // Inputs:
+    //   V  #V by 3 list of vertex positions
+    //   F  #F by 3 list of triangle indices into V
+    //   I  #I list of facets to consider
+    // Outputs:
+    //   v_index  index of outer vertex
+    //   A  #A list of facets incident to the outer vertex
+    template <
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename IndexType,
+        typename DerivedA
+        >
+    IGL_INLINE void outer_vertex(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::PlainObjectBase<DerivedF> & F,
+            const Eigen::PlainObjectBase<DerivedI> & I,
+            IndexType & v_index,
+            Eigen::PlainObjectBase<DerivedA> & A);
+    // Find an edge that is reachable from infinity without crossing any faces.
+    // Such edge is called "outer edge."
+    //
+    // Precondition: The input mesh must have all self-intersection resolved
+    // and no duplicated vertices.  The correctness of the output depends on
+    // the fact that there is no edge overlap.  See
+    // cgal::remesh_self_intersections.h for how to obtain such input.
+    //
+    // Inputs:
+    //   V  #V by 3 list of vertex positions
+    //   F  #F by 3 list of triangle indices into V
+    //   I  #I list of facets to consider
+    // Outputs:
+    //   v1 index of the first end point of outer edge
+    //   v2 index of the second end point of outer edge
+    //   A  #A list of facets incident to the outer edge
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename IndexType,
+        typename DerivedA
+        >
+    IGL_INLINE void outer_edge(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::PlainObjectBase<DerivedF> & F,
+            const Eigen::PlainObjectBase<DerivedI> & I,
+            IndexType & v1,
+            IndexType & v2,
+            Eigen::PlainObjectBase<DerivedA> & A);
+
+
+    // Find a facet that is reachable from infinity without crossing any faces.
+    // Such facet is called "outer facet."
+    //
+    // Precondition: The input mesh must have all self-intersection resolved.
+    // I.e there is no duplicated vertices, no overlapping edge and no
+    // intersecting faces (the only exception is there could be topologically
+    // duplicated faces).  See cgal::remesh_self_intersections.h for how to
+    // obtain such input.
+    //
+    // Inputs:
+    //   V  #V by 3 list of vertex positions
+    //   F  #F by 3 list of triangle indices into V
+    //   N  #N by 3 list of face normals
+    //   I  #I list of facets to consider
+    // Outputs:
+    //   f  Index of the outer facet.
+    //   flipped  true iff the normal of f points inwards.
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedN,
+        typename DerivedI,
+        typename IndexType
+        >
+    IGL_INLINE void outer_facet(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::PlainObjectBase<DerivedF> & F,
+            const Eigen::PlainObjectBase<DerivedN> & N,
+            const Eigen::PlainObjectBase<DerivedI> & I,
+            IndexType & f,
+            bool & flipped);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "outer_element.cpp"
+#endif
+#endif

include/igl/cgal/outer_facet.cpp

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "outer_facet.h"
-#include "../outer_element.h"
+#include "outer_element.h"
 #include "order_facets_around_edge.h"
 #include <algorithm>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
@@ -86,4 +86,6 @@ IGL_INLINE void igl::cgal::outer_facet(
 template void igl::cgal::outer_facet<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>, int>(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&, int&, bool&);
 template void igl::cgal::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, unsigned long&, bool&);
+template void igl::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
 #endif

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,261 @@
+// 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);
+    //}
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/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.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/cgal/remesh_self_intersections.cpp

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

include/igl/cgal/signed_distance_isosurface.cpp

@@ -16,7 +16,6 @@
 #include "../centroid.h"
 #include "../WindingNumberAABB.h"
 #include "../matlab_format.h"
-#include "../remove_unreferenced.h"
 
 #include <CGAL/Surface_mesh_default_triangulation_3.h>
 #include <CGAL/Complex_2_in_triangulation_3.h>

include/igl/copyleft/marching_cubes.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_COPYLEFT_MARCHINGCUBES_H
 #define IGL_COPYLEFT_MARCHINGCUBES_H
-#include "igl_inline.h"
+#include "../igl_inline.h"
 
 #include <Eigen/Core>
 namespace igl 

include/igl/extract_manifold_patches.cpp

@@ -0,0 +1,72 @@
+#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;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template unsigned long igl::extract_manifold_patches<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

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

include/igl/facet_components.cpp

@@ -87,4 +87,5 @@ IGL_INLINE void igl::facet_components(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::facet_components<long, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::facet_components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::__1::vector<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >, std::__1::allocator<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/remove_unreferenced.cpp

@@ -102,4 +102,5 @@ IGL_INLINE void igl::remove_unreferenced(
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, 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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > 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> >&);
 #endif

include/igl/triangle_triangle_adjacency.cpp

@@ -241,4 +241,5 @@ template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0,
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::triangle_triangle_adjacency<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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, bool, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, int, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > > const&, bool, std::__1::vector<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >, std::__1::allocator<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > > > >&, std::__1::vector<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >, std::__1::allocator<std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > > > >&);
 #endif

include/igl/unique_edge_map.cpp

@@ -51,4 +51,5 @@ template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen:
 template void igl::unique_edge_map<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::Matrix<int, -1, 1, 0, -1, 1>, int>(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<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > >&);
 #endif 

tutorial/609_Boolean/main.cpp

@@ -26,7 +26,7 @@ const char * MESH_BOOLEAN_TYPE_NAMES[] =
 
 void update(igl::viewer::Viewer &viewer)
 {
-  igl::boolean::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC,J,I);
+  igl::boolean::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC,J);
   Eigen::MatrixXd C(FC.rows(),3);
   for(size_t f = 0;f<C.rows();f++)
   {

tutorial/705_MarchingCubes/main.cpp

@@ -1,4 +1,4 @@
-#include <igl/marching_cubes.h>
+#include <igl/copyleft/marching_cubes.h>
 #include <igl/signed_distance.h>
 #include <igl/read_triangle_mesh.h>
 #include <igl/viewer/Viewer.h>
@@ -55,8 +55,8 @@ int main(int argc, char * argv[])
   cout<<"Marching cubes..."<<endl;
   MatrixXd SV,BV;
   MatrixXi SF,BF;
-  marching_cubes(S,GV,res(0),res(1),res(2),SV,SF);
-  marching_cubes(B,GV,res(0),res(1),res(2),BV,BF);
+  igl::copyleft::marching_cubes(S,GV,res(0),res(1),res(2),SV,SF);
+  igl::copyleft::marching_cubes(B,GV,res(0),res(1),res(2),BV,BF);
 
   cout<<R"(Usage:
 '1'  Show original mesh.