merge

Former-commit-id: aaad45044a53687f1cd502d91d3a619e414ad2f2

include/igl/colon.cpp

@@ -52,4 +52,5 @@ template void igl::colon<int, long, int>(int, long, Eigen::Matrix<int, -1, 1, 0,
 template void igl::colon<int, int, int>(int, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 template void igl::colon<int,long long int,int>(int,long long int,Eigen::Matrix<int,-1,1,0,-1,1> &);
 template void igl::colon<int, int, int, int>(int, int, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
 #endif

include/igl/copyleft/boolean/BinaryWindingNumberOperations.h

@@ -20,110 +20,110 @@ namespace igl
   {
     namespace boolean
     {
-        template <igl::copyleft::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 <igl::copyleft::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::copyleft::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<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::copyleft::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<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::copyleft::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<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::copyleft::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<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::copyleft::boolean::MESH_BOOLEAN_TYPE_RESOLVE> {
-            public:
-                template<typename DerivedW>
-                    typename DerivedW::Scalar operator()(
-                            const Eigen::PlainObjectBase<DerivedW>& /*win_nums*/) const {
-                        return true;
-                    }
-        };
+      template <>
+      class BinaryWindingNumberOperations<igl::copyleft::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;
+      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
-        };
+      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<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_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;
-                    }
-        };
+      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;
+      typedef WindingNumberFilter<KEEP_INSIDE> KeepInside;
+      typedef WindingNumberFilter<KEEP_ALL> KeepAll;
     }
   }
 }

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -15,274 +15,133 @@
 #include "../../remove_unreferenced.h"
 #include "../../unique_simplices.h"
 #include "../../slice.h"
+#include "../../resolve_duplicated_faces.h"
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <algorithm>
 
-
-namespace igl {
-  namespace copyleft {
-    namespace boolean {
-        namespace mesh_boolean_helper {
-            typedef CGAL::Epeck Kernel;
-            typedef Kernel::FT ExactScalar;
-
-            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::copyleft::cgal::RemeshSelfIntersectionsParam params;
-
-                DerivedVo Vr;
-                DerivedFo Fr;
-                Eigen::MatrixXi IF;
-                igl::copyleft::cgal::remesh_self_intersections(V, F, params, Vr, Fr, IF, J, I);
-                assert(I.size() == Vr.rows());
-
-                // Merge coinciding vertices into non-manifold vertices.
-                std::for_each(Fr.data(), Fr.data()+Fr.size(),
-                        [&I](typename DerivedF::Scalar& a) { a=I[a]; });
-
-                // 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((size_t) 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]);
-                }
-            }
-        }
-    }
-  }
-}
-
 template <
-typename DerivedVA,
-typename DerivedFA,
-typename DerivedVB,
-typename DerivedFB,
-typename WindingNumberOp,
-typename KeepFunc,
-typename ResolveFunc,
-typename DerivedVC,
-typename DerivedFC,
-typename DerivedJ>
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename WindingNumberOp,
+  typename KeepFunc,
+  typename ResolveFunc,
+  typename DerivedVC,
+  typename DerivedFC,
+  typename DerivedJ>
 IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA> & VA,
-  const Eigen::PlainObjectBase<DerivedFA> & FA,
-  const Eigen::PlainObjectBase<DerivedVB> & VB,
-  const Eigen::PlainObjectBase<DerivedFB> & FB,
-  const 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::copyleft::boolean::mesh_boolean_helper;
+    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) {
 
   typedef typename DerivedVC::Scalar Scalar;
   //typedef typename DerivedFC::Scalar Index;
+  typedef CGAL::Epeck Kernel;
+  typedef Kernel::FT ExactScalar;
   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;
+    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);
+  {
+      DerivedVA VV(VA.rows() + VB.rows(), 3);
+      DerivedFC FF(FA.rows() + FB.rows(), 3);
+      VV << VA, VB;
+      FF << FA, FB.array() + VA.rows();
+      resolve_fun(VV, FF, 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; });
+      [&](int i) { return i<FA.rows() ? 0:1; });
   igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
   assert((size_t)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;
+    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);
+    assert(W.cols() == 4);
   }
 
   // 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);
+    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);
   }
 
   // Extract boundary separating inside from outside.
   auto index_to_signed_index = [&](size_t i, bool ori) -> int{
-      return (i+1)*(ori?1:-1);
+    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));
-      }
+    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);
+  DerivedJ  kept_face_indices(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];
+    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];
   }
 
-
   // Finally, remove duplicated faces and unreferenced vertices.
   {
-      DerivedFC G;
-      DerivedJ J;
-      resolve_duplicated_faces(kept_faces, kept_face_indices, G, J);
+    DerivedFC G;
+    DerivedJ JJ;
+    igl::resolve_duplicated_faces(kept_faces, G, JJ);
+    igl::slice(kept_face_indices, JJ, 1, J);
 
-      MatrixX3S Vs(V.rows(), V.cols());
-      for (size_t i=0; i<(size_t)V.rows(); i++) {
-          for (size_t j=0; j<(size_t)V.cols(); j++) {
-              igl::copyleft::cgal::assign_scalar(V(i,j), Vs(i,j));
-          }
+    MatrixX3S Vs(V.rows(), V.cols());
+    for (size_t i=0; i<(size_t)V.rows(); i++)
+    {
+      for (size_t j=0; j<(size_t)V.cols(); j++)
+      {
+        igl::copyleft::cgal::assign_scalar(V(i,j), Vs(i,j));
       }
-      Eigen::VectorXi newIM;
-      igl::remove_unreferenced(Vs,G,VC,FC,newIM);
+    }
+    Eigen::VectorXi newIM;
+    igl::remove_unreferenced(Vs,G,VC,FC,newIM);
   }
 }
 
@@ -296,48 +155,54 @@ template <
   typename DerivedFC,
   typename DerivedJ>
 IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
-  const Eigen::PlainObjectBase<DerivedVA > & VA,
-  const Eigen::PlainObjectBase<DerivedFA > & FA,
-  const Eigen::PlainObjectBase<DerivedVB > & VB,
-  const Eigen::PlainObjectBase<DerivedFB > & FB,
-  const MeshBooleanType & type,
-  const ResolveFunc& resolve_func,
-  Eigen::PlainObjectBase<DerivedVC > & VC,
-  Eigen::PlainObjectBase<DerivedFC > & FC,
-  Eigen::PlainObjectBase<DerivedJ > & J)
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    const ResolveFunc& resolve_func,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC,
+    Eigen::PlainObjectBase<DerivedJ > & J)
 {
-  using namespace igl::copyleft::boolean::mesh_boolean_helper;
+  typedef CGAL::Epeck Kernel;
+  typedef Kernel::FT ExactScalar;
+  typedef Eigen::Matrix<
+    ExactScalar,
+    Eigen::Dynamic,
+    Eigen::Dynamic,
+    DerivedVC::IsRowMajor> MatrixXES;
 
-    switch (type) {
-        case MESH_BOOLEAN_TYPE_UNION:
-            igl::copyleft::boolean::mesh_boolean(
-                    VA, FA, VB, FB, igl::copyleft::boolean::BinaryUnion(),
-                    igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-            break;
-        case MESH_BOOLEAN_TYPE_INTERSECT:
-            igl::copyleft::boolean::mesh_boolean(
-                    VA, FA, VB, FB, igl::copyleft::boolean::BinaryIntersect(),
-                    igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-            break;
-        case MESH_BOOLEAN_TYPE_MINUS:
-            igl::copyleft::boolean::mesh_boolean(
-                    VA, FA, VB, FB, igl::copyleft::boolean::BinaryMinus(),
-                    igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-            break;
-        case MESH_BOOLEAN_TYPE_XOR:
-            igl::copyleft::boolean::mesh_boolean(
-                    VA, FA, VB, FB, igl::copyleft::boolean::BinaryXor(),
-                    igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
-            break;
-        case MESH_BOOLEAN_TYPE_RESOLVE:
-            //op = binary_resolve();
-            igl::copyleft::boolean::mesh_boolean(
-                    VA, FA, VB, FB, igl::copyleft::boolean::BinaryResolve(),
-                    igl::copyleft::boolean::KeepAll(), resolve_func, VC, FC, J);
-            break;
-        default:
-            throw std::runtime_error("Unsupported boolean type.");
-    }
+  switch (type) {
+    case MESH_BOOLEAN_TYPE_UNION:
+      igl::copyleft::boolean::mesh_boolean(
+          VA, FA, VB, FB, igl::copyleft::boolean::BinaryUnion(),
+          igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
+      break;
+    case MESH_BOOLEAN_TYPE_INTERSECT:
+      igl::copyleft::boolean::mesh_boolean(
+          VA, FA, VB, FB, igl::copyleft::boolean::BinaryIntersect(),
+          igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
+      break;
+    case MESH_BOOLEAN_TYPE_MINUS:
+      igl::copyleft::boolean::mesh_boolean(
+          VA, FA, VB, FB, igl::copyleft::boolean::BinaryMinus(),
+          igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
+      break;
+    case MESH_BOOLEAN_TYPE_XOR:
+      igl::copyleft::boolean::mesh_boolean(
+          VA, FA, VB, FB, igl::copyleft::boolean::BinaryXor(),
+          igl::copyleft::boolean::KeepInside(), resolve_func, VC, FC, J);
+      break;
+    case MESH_BOOLEAN_TYPE_RESOLVE:
+      //op = binary_resolve();
+      igl::copyleft::boolean::mesh_boolean(
+          VA, FA, VB, FB, igl::copyleft::boolean::BinaryResolve(),
+          igl::copyleft::boolean::KeepAll(), resolve_func, VC, FC, J);
+      break;
+    default:
+      throw std::runtime_error("Unsupported boolean type.");
+  }
 }
 
 template <
@@ -358,19 +223,44 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
   Eigen::PlainObjectBase<DerivedFC > & FC,
   Eigen::PlainObjectBase<DerivedJ > & J)
 {
-  using namespace igl::copyleft::boolean::mesh_boolean_helper;
+  typedef CGAL::Epeck Kernel;
+  typedef Kernel::FT ExactScalar;
   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>;
+      const Eigen::PlainObjectBase<DerivedVA>&,
+      const Eigen::PlainObjectBase<DerivedFA>&,
+      Eigen::PlainObjectBase<MatrixXES>&,
+      Eigen::PlainObjectBase<DerivedFC>&,
+      Eigen::PlainObjectBase<DerivedJ>&)> resolve_func =
+    [](const Eigen::PlainObjectBase<DerivedVA>& V,
+        const Eigen::PlainObjectBase<DerivedFA>& F,
+        Eigen::PlainObjectBase<MatrixXES>& Vo,
+        Eigen::PlainObjectBase<DerivedFC>& Fo,
+        Eigen::PlainObjectBase<DerivedJ>& J) {
+      Eigen::VectorXi I;
+      igl::copyleft::cgal::RemeshSelfIntersectionsParam params;
+
+      MatrixXES Vr;
+      DerivedFC Fr;
+      Eigen::MatrixXi IF;
+      igl::copyleft::cgal::remesh_self_intersections(
+          V, F, params, Vr, Fr, IF, J, I);
+      assert(I.size() == Vr.rows());
+
+      // Merge coinciding vertices into non-manifold vertices.
+      std::for_each(Fr.data(), Fr.data()+Fr.size(),
+          [&I](typename DerivedFC::Scalar& a) { a=I[a]; });
+
+      // Remove unreferenced vertices.
+      Eigen::VectorXi UIM;
+      igl::remove_unreferenced(Vr, Fr, Vo, Fo, UIM);
+    };
+
   return mesh_boolean(VA,FA,VB,FB,type,resolve_func,VC,FC,J);
 }
 
@@ -382,15 +272,15 @@ typename DerivedFB,
 typename DerivedVC,
 typename DerivedFC>
 IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
-        const Eigen::PlainObjectBase<DerivedVA > & VA,
-        const Eigen::PlainObjectBase<DerivedFA > & FA,
-        const Eigen::PlainObjectBase<DerivedVB > & VB,
-        const Eigen::PlainObjectBase<DerivedFB > & FB,
-        const MeshBooleanType & type,
-        Eigen::PlainObjectBase<DerivedVC > & VC,
-        Eigen::PlainObjectBase<DerivedFC > & FC) {
-    Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
-    return igl::copyleft::boolean::mesh_boolean(VA,FA,VB,FB,type,VC,FC,J);
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC) {
+  Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
+  return igl::copyleft::boolean::mesh_boolean(VA,FA,VB,FB,type,VC,FC,J);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/copyleft/boolean/mesh_boolean.h

@@ -23,7 +23,7 @@ namespace igl
     {
       //  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
@@ -42,7 +42,6 @@ namespace igl
       //
       //  See also: mesh_boolean_cork, intersect_other,
       //  remesh_self_intersections
-      //     
       template <
         typename DerivedVA,
         typename DerivedFA,
@@ -55,16 +54,17 @@ namespace igl
         typename DerivedFC,
         typename DerivedJ>
       IGL_INLINE void mesh_boolean(
-        const Eigen::PlainObjectBase<DerivedVA> & VA,
-        const Eigen::PlainObjectBase<DerivedFA> & FA,
-        const Eigen::PlainObjectBase<DerivedVB> & VB,
-        const Eigen::PlainObjectBase<DerivedFB> & FB,
-        const WindingNumberOp& wind_num_op,
-        const KeepFunc& keep,
-        const ResolveFunc& resolve_fun,
-        Eigen::PlainObjectBase<DerivedVC > & VC,
-        Eigen::PlainObjectBase<DerivedFC > & FC,
-        Eigen::PlainObjectBase<DerivedJ > & J);
+          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);
+
       //  Inputs:
       //    VA  #VA by 3 list of vertex positions of first mesh
       //    FA  #FA by 3 list of triangle indices into VA
@@ -90,15 +90,16 @@ namespace igl
         typename DerivedFC,
         typename DerivedJ>
       IGL_INLINE void mesh_boolean(
-        const Eigen::PlainObjectBase<DerivedVA > & VA,
-        const Eigen::PlainObjectBase<DerivedFA > & FA,
-        const Eigen::PlainObjectBase<DerivedVB > & VB,
-        const Eigen::PlainObjectBase<DerivedFB > & FB,
-        const MeshBooleanType & type,
+          const Eigen::PlainObjectBase<DerivedVA > & VA,
+          const Eigen::PlainObjectBase<DerivedFA > & FA,
+          const Eigen::PlainObjectBase<DerivedVB > & VB,
+          const Eigen::PlainObjectBase<DerivedFB > & FB,
+          const MeshBooleanType & type,
         const ResolveFunc& resolve_func,
-        Eigen::PlainObjectBase<DerivedVC > & VC,
-        Eigen::PlainObjectBase<DerivedFC > & FC,
-        Eigen::PlainObjectBase<DerivedJ > & J);
+          Eigen::PlainObjectBase<DerivedVC > & VC,
+          Eigen::PlainObjectBase<DerivedFC > & FC,
+          Eigen::PlainObjectBase<DerivedJ > & J);
+
       //  Inputs:
       //    VA  #VA by 3 list of vertex positions of first mesh
       //    FA  #FA by 3 list of triangle indices into VA
@@ -129,6 +130,7 @@ namespace igl
         Eigen::PlainObjectBase<DerivedVC > & VC,
         Eigen::PlainObjectBase<DerivedFC > & FC,
         Eigen::PlainObjectBase<DerivedJ > & J);
+
       //  Inputs:
       //    VA  #VA by 3 list of vertex positions of first mesh
       //    FA  #FA by 3 list of triangle indices into VA
@@ -146,13 +148,13 @@ namespace igl
         typename DerivedVC,
         typename DerivedFC>
       IGL_INLINE void mesh_boolean(
-        const Eigen::PlainObjectBase<DerivedVA > & VA,
-        const Eigen::PlainObjectBase<DerivedFA > & FA,
-        const Eigen::PlainObjectBase<DerivedVB > & VB,
-        const Eigen::PlainObjectBase<DerivedFB > & FB,
-        const MeshBooleanType & type,
-        Eigen::PlainObjectBase<DerivedVC > & VC,
-        Eigen::PlainObjectBase<DerivedFC > & FC);
+          const Eigen::PlainObjectBase<DerivedVA > & VA,
+          const Eigen::PlainObjectBase<DerivedFA > & FA,
+          const Eigen::PlainObjectBase<DerivedVB > & VB,
+          const Eigen::PlainObjectBase<DerivedFB > & FB,
+          const MeshBooleanType & type,
+          Eigen::PlainObjectBase<DerivedVC > & VC,
+          Eigen::PlainObjectBase<DerivedFC > & FC);
     }
   }
 }

include/igl/copyleft/cgal/order_facets_around_edge.cpp

@@ -3,40 +3,11 @@
 
 #include <stdexcept>
 
-namespace igl
-{
-  namespace copyleft
-  {
-    namespace cgal
-    {
-        namespace order_facets_around_edges_helper
-        {
-            template<typename T>
-            std::vector<size_t> index_sort(const std::vector<T>& data)
-            {
-                const size_t len = data.size();
-                std::vector<size_t> order(len);
-                for (size_t i=0; i<len; i++) 
-                {
-                    order[i] = i;
-                }
-                auto comp = [&](size_t i, size_t j) 
-                {
-                    return data[i] < data[j];
-                };
-                std::sort(order.begin(), order.end(), comp);
-                return order;
-            }
-        }
-    }
-  }
-}
-
 // adj_faces contains signed index starting from +- 1.
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedI >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI >
 void igl::copyleft::cgal::order_facets_around_edge(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
@@ -45,254 +16,261 @@ void igl::copyleft::cgal::order_facets_around_edge(
     const std::vector<int>& adj_faces,
     Eigen::PlainObjectBase<DerivedI>& order, bool debug)
 {
-    using namespace igl::copyleft::cgal::order_facets_around_edges_helper;
-
-    // Although we only need exact predicates in the algorithm,
-    // exact constructions are needed to avoid degeneracies due to
-    // casting to double.
-    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
-    typedef K::Point_3 Point_3;
-    typedef K::Plane_3 Plane_3;
-
-    auto get_face_index = [&](int adj_f)->size_t
-    {
-        return abs(adj_f) - 1;
-    };
-
-    auto get_opposite_vertex = [&](size_t fid)->size_t
-    {
-        typedef typename DerivedF::Scalar Index;
-        if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
-        if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
-        if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
-        assert(false);
-        return -1;
-    };
+  // Although we only need exact predicates in the algorithm,
+  // exact constructions are needed to avoid degeneracies due to
+  // casting to double.
+  typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+  typedef K::Point_3 Point_3;
+  typedef K::Plane_3 Plane_3;
+
+  auto get_face_index = [&](int adj_f)->size_t
+  {
+    return abs(adj_f) - 1;
+  };
 
-    // Handle base cases
-    if (adj_faces.size() == 0) 
-    {
-        order.resize(0, 1);
-        return;
-    } else if (adj_faces.size() == 1)
+  auto get_opposite_vertex = [&](size_t fid)->size_t
+  {
+    typedef typename DerivedF::Scalar Index;
+    if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
+    if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
+    if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
+    assert(false);
+    return -1;
+  };
+
+  // Handle base cases
+  if (adj_faces.size() == 0) 
+  {
+    order.resize(0, 1);
+    return;
+  } else if (adj_faces.size() == 1)
+  {
+    order.resize(1, 1);
+    order(0, 0) = 0;
+    return;
+  } else if (adj_faces.size() == 2)
+  {
+    const size_t o1 = get_opposite_vertex(get_face_index(adj_faces[0]));
+    const size_t o2 = get_opposite_vertex(get_face_index(adj_faces[1]));
+    const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
+    const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
+    const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
+    const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
+    order.resize(2, 1);
+    switch (CGAL::orientation(ps, pd, p1, p2))
     {
-        order.resize(1, 1);
+      case CGAL::POSITIVE:
+        order(0, 0) = 1;
+        order(1, 0) = 0;
+        break;
+      case CGAL::NEGATIVE:
         order(0, 0) = 0;
-        return;
-    } else if (adj_faces.size() == 2)
-    {
-        const size_t o1 = get_opposite_vertex(get_face_index(adj_faces[0]));
-        const size_t o2 = get_opposite_vertex(get_face_index(adj_faces[1]));
-        const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
-        const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
-        const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
-        const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
-        order.resize(2, 1);
-        switch (CGAL::orientation(ps, pd, p1, p2))
+        order(1, 0) = 1;
+        break;
+      case CGAL::COPLANAR:
         {
+          switch (CGAL::coplanar_orientation(ps, pd, p1, p2)) {
             case CGAL::POSITIVE:
-                order(0, 0) = 1;
-                order(1, 0) = 0;
-                break;
+              // Duplicated face, use index to break tie.
+              order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
+              order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+              break;
             case CGAL::NEGATIVE:
-                order(0, 0) = 0;
-                order(1, 0) = 1;
-                break;
-            case CGAL::COPLANAR:
-                {
-                    switch (CGAL::coplanar_orientation(ps, pd, p1, p2)) {
-                        case CGAL::POSITIVE:
-                            // Duplicated face, use index to break tie.
-                            order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
-                            order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
-                            break;
-                        case CGAL::NEGATIVE:
-                            // Coplanar faces, one on each side of the edge.
-                            // It is equally valid to order them (0, 1) or (1, 0).
-                            // I cannot think of any reason to prefer one to the
-                            // other.  So just use (0, 1) ordering by default.
-                            order(0, 0) = 0;
-                            order(1, 0) = 1;
-                            break;
-                        case CGAL::COLLINEAR:
-                            std::cerr << "Degenerated triangle detected." <<
-                                std::endl;
-                            assert(false);
-                            break;
-                        default:
-                            assert(false);
-                    }
-                }
-                break;
+              // Coplanar faces, one on each side of the edge.
+              // It is equally valid to order them (0, 1) or (1, 0).
+              // I cannot think of any reason to prefer one to the
+              // other.  So just use (0, 1) ordering by default.
+              order(0, 0) = 0;
+              order(1, 0) = 1;
+              break;
+            case CGAL::COLLINEAR:
+              std::cerr << "Degenerated triangle detected." <<
+                std::endl;
+              assert(false);
+              break;
             default:
-                assert(false);
+              assert(false);
+          }
         }
-        return;
+        break;
+      default:
+        assert(false);
     }
+    return;
+  }
 
-    const size_t num_adj_faces = adj_faces.size();
-    const size_t o = get_opposite_vertex( get_face_index(adj_faces[0]));
-    const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
-    const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
-    const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
-    const Plane_3 separator(p_s, p_d, p_o);
-    if (separator.is_degenerate()) {
-        throw std::runtime_error(
-                "Cannot order facets around edge due to degenerated facets");
-    }
+  const size_t num_adj_faces = adj_faces.size();
+  const size_t o = get_opposite_vertex( get_face_index(adj_faces[0]));
+  const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
+  const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
+  const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
+  const Plane_3 separator(p_s, p_d, p_o);
+  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++)
-    {
-        const size_t o = get_opposite_vertex( get_face_index(adj_faces[i]));
-        opposite_vertices.emplace_back(
-                V(o, 0), V(o, 1), V(o, 2));
-    }
+  std::vector<Point_3> opposite_vertices;
+  for (size_t i=0; i<num_adj_faces; i++)
+  {
+    const size_t o = get_opposite_vertex( get_face_index(adj_faces[i]));
+    opposite_vertices.emplace_back(
+        V(o, 0), V(o, 1), V(o, 2));
+  }
 
-    std::vector<int> positive_side;
-    std::vector<int> negative_side;
-    std::vector<int> tie_positive_oriented;
-    std::vector<int> tie_negative_oriented;
+  std::vector<int> positive_side;
+  std::vector<int> negative_side;
+  std::vector<int> tie_positive_oriented;
+  std::vector<int> tie_negative_oriented;
 
-    std::vector<size_t> positive_side_index;
-    std::vector<size_t> negative_side_index;
-    std::vector<size_t> tie_positive_oriented_index;
-    std::vector<size_t> tie_negative_oriented_index;
+  std::vector<size_t> positive_side_index;
+  std::vector<size_t> negative_side_index;
+  std::vector<size_t> tie_positive_oriented_index;
+  std::vector<size_t> tie_negative_oriented_index;
 
-    for (size_t i=0; i<num_adj_faces; i++)
-    {
-        const int f = adj_faces[i];
-        const Point_3& p_a = opposite_vertices[i];
-        auto orientation = separator.oriented_side(p_a);
-        switch (orientation) {
-            case CGAL::ON_POSITIVE_SIDE:
-                positive_side.push_back(f);
-                positive_side_index.push_back(i);
-                break;
-            case CGAL::ON_NEGATIVE_SIDE:
-                negative_side.push_back(f);
-                negative_side_index.push_back(i);
-                break;
-            case CGAL::ON_ORIENTED_BOUNDARY:
-                {
-                    auto inplane_orientation = CGAL::coplanar_orientation(
-                            p_s, p_d, p_o, p_a);
-                    switch (inplane_orientation) {
-                        case CGAL::POSITIVE:
-                            tie_positive_oriented.push_back(f);
-                            tie_positive_oriented_index.push_back(i);
-                            break;
-                        case CGAL::NEGATIVE:
-                            tie_negative_oriented.push_back(f);
-                            tie_negative_oriented_index.push_back(i);
-                            break;
-                        case CGAL::COLLINEAR:
-                        default:
-                            throw std::runtime_error(
-                                    "Degenerated facet detected.");
-                            break;
-                    }
-                }
-                break;
+  for (size_t i=0; i<num_adj_faces; i++)
+  {
+    const int f = adj_faces[i];
+    const Point_3& p_a = opposite_vertices[i];
+    auto orientation = separator.oriented_side(p_a);
+    switch (orientation) {
+      case CGAL::ON_POSITIVE_SIDE:
+        positive_side.push_back(f);
+        positive_side_index.push_back(i);
+        break;
+      case CGAL::ON_NEGATIVE_SIDE:
+        negative_side.push_back(f);
+        negative_side_index.push_back(i);
+        break;
+      case CGAL::ON_ORIENTED_BOUNDARY:
+        {
+          auto inplane_orientation = CGAL::coplanar_orientation(
+              p_s, p_d, p_o, p_a);
+          switch (inplane_orientation) {
+            case CGAL::POSITIVE:
+              tie_positive_oriented.push_back(f);
+              tie_positive_oriented_index.push_back(i);
+              break;
+            case CGAL::NEGATIVE:
+              tie_negative_oriented.push_back(f);
+              tie_negative_oriented_index.push_back(i);
+              break;
+            case CGAL::COLLINEAR:
             default:
-                // Should not be here.
-                throw std::runtime_error("Unknown CGAL state detected.");
+              throw std::runtime_error(
+                  "Degenerated facet detected.");
+              break;
+          }
         }
+        break;
+      default:
+        // Should not be here.
+        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;
+  }
+  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, 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);
-
-    // Copy results into order vector.
-    const size_t tie_positive_size = tie_positive_oriented.size();
-    const size_t tie_negative_size = tie_negative_oriented.size();
-    const size_t positive_size = positive_order.size();
-    const size_t negative_size = negative_order.size();
-
-    order.resize(
+  auto index_sort = [](std::vector<int>& data) -> std::vector<size_t>{
+    const size_t len = data.size();
+    std::vector<size_t> order(len);
+    for (size_t i=0; i<len; i++) { order[i] = i; }
+    auto comp = [&](size_t i, size_t j) { return data[i] < data[j]; };
+    std::sort(order.begin(), order.end(), comp);
+    return order;
+  };
+
+  Eigen::PlainObjectBase<DerivedI> positive_order, 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);
+
+  // Copy results into order vector.
+  const size_t tie_positive_size = tie_positive_oriented.size();
+  const size_t tie_negative_size = tie_negative_oriented.size();
+  const size_t positive_size = positive_order.size();
+  const size_t negative_size = negative_order.size();
+
+  order.resize(
       tie_positive_size + positive_size + tie_negative_size + negative_size,1);
 
-    size_t count=0;
-    for (size_t i=0; i<tie_positive_size; i++)
-    {
-        order(count+i, 0) = tie_positive_oriented_index[tie_positive_order[i]];
-    }
-    count += tie_positive_size;
+  size_t count=0;
+  for (size_t i=0; i<tie_positive_size; i++)
+  {
+    order(count+i, 0) = tie_positive_oriented_index[tie_positive_order[i]];
+  }
+  count += tie_positive_size;
 
-    for (size_t i=0; i<negative_size; i++) 
-    {
-        order(count+i, 0) = negative_side_index[negative_order(i, 0)];
-    }
-    count += negative_size;
+  for (size_t i=0; i<negative_size; i++) 
+  {
+    order(count+i, 0) = negative_side_index[negative_order(i, 0)];
+  }
+  count += negative_size;
 
-    for (size_t i=0; i<tie_negative_size; i++)
-    {
-        order(count+i, 0) = tie_negative_oriented_index[tie_negative_order[i]];
-    }
-    count += tie_negative_size;
+  for (size_t i=0; i<tie_negative_size; i++)
+  {
+    order(count+i, 0) = tie_negative_oriented_index[tie_negative_order[i]];
+  }
+  count += tie_negative_size;
 
-    for (size_t i=0; i<positive_size; i++)
-    {
-        order(count+i, 0) = positive_side_index[positive_order(i, 0)];
-    }
-    count += positive_size;
-    assert(count == num_adj_faces);
+  for (size_t i=0; i<positive_size; i++)
+  {
+    order(count+i, 0) = positive_side_index[positive_order(i, 0)];
+  }
+  count += positive_size;
+  assert(count == num_adj_faces);
 
-    // Find the correct start point.
-    size_t start_idx = 0;
-    for (size_t i=0; i<num_adj_faces; i++)
+  // Find the correct start point.
+  size_t start_idx = 0;
+  for (size_t i=0; i<num_adj_faces; i++)
+  {
+    const Point_3& p_a = opposite_vertices[order(i, 0)];
+    const Point_3& p_b =
+      opposite_vertices[order((i+1)%num_adj_faces, 0)];
+    auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
+    if (orientation == CGAL::POSITIVE)
     {
-        const Point_3& p_a = opposite_vertices[order(i, 0)];
-        const Point_3& p_b =
-            opposite_vertices[order((i+1)%num_adj_faces, 0)];
-        auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
-        if (orientation == CGAL::POSITIVE)
-        {
-            // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
-            // more than 180 degrees.
-            start_idx = (i+1)%num_adj_faces;
-            break;
-        } else if (orientation == CGAL::COPLANAR &&
-                Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
-                Plane_3(p_s, p_d, p_b).orthogonal_direction())
-        {
-            // All 4 points are coplanar, but p_a and p_b are on each side of
-            // the edge (p_s, p_d).  This means the angle between triangle
-            // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
-            start_idx = (i+1)%num_adj_faces;
-            break;
-        }
-    }
-    DerivedI circular_order = order;
-    for (size_t i=0; i<num_adj_faces; i++)
+      // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
+      // more than 180 degrees.
+      start_idx = (i+1)%num_adj_faces;
+      break;
+    } else if (orientation == CGAL::COPLANAR &&
+        Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
+        Plane_3(p_s, p_d, p_b).orthogonal_direction())
     {
-        order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
+      // All 4 points are coplanar, but p_a and p_b are on each side of
+      // the edge (p_s, p_d).  This means the angle between triangle
+      // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
+      start_idx = (i+1)%num_adj_faces;
+      break;
     }
+  }
+  DerivedI circular_order = order;
+  for (size_t i=0; i<num_adj_faces; i++)
+  {
+    order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
+  }
 }
 
 template<

include/igl/copyleft/cgal/order_facets_around_edge.h

@@ -16,54 +16,54 @@ namespace igl {
   namespace copyleft
   {
     namespace cgal {
-        // Given a directed edge, sort its adjacent faces.  Assuming the
-        // directed edge is (s, d).  Sort the adjacent faces clockwise around the
-        // axis (d - s), i.e. left-hand rule.  An adjacent face is consistently
-        // oriented if it contains (d, s) as a directed edge.
-        //
-        // For overlapping faces, break the tie using signed face index, smaller
-        // signed index comes before the larger signed index.  Signed index is
-        // computed as (consistent? 1:-1) * (face_index + 1).
-        //
-        // Inputs:
-        //   V          #V by 3 list of vertices.
-        //   F          #F by 3 list of faces
-        //   s          Index of source vertex.
-        //   d          Index of desination vertex.
-        //   adj_faces  List of adjacent face signed indices.
-        //
-        // Output:
-        //   order      List of face indices that orders adjacent faces around
-        //              edge (s, d) clockwise.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedI >
-        IGL_INLINE
-        void order_facets_around_edge(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                size_t s, size_t d, 
-                const std::vector<int>& adj_faces,
-                Eigen::PlainObjectBase<DerivedI>& order,
-                bool debug=false);
+      // Given a directed edge, sort its adjacent faces.  Assuming the
+      // directed edge is (s, d).  Sort the adjacent faces clockwise around the
+      // axis (d - s), i.e. left-hand rule.  An adjacent face is consistently
+      // oriented if it contains (d, s) as a directed edge.
+      //
+      // For overlapping faces, break the tie using signed face index, smaller
+      // signed index comes before the larger signed index.  Signed index is
+      // computed as (consistent? 1:-1) * (face_index + 1).
+      //
+      // Inputs:
+      //   V          #V by 3 list of vertices.
+      //   F          #F by 3 list of faces
+      //   s          Index of source vertex.
+      //   d          Index of desination vertex.
+      //   adj_faces  List of adjacent face signed indices.
+      //
+      // Output:
+      //   order      List of face indices that orders adjacent faces around
+      //              edge (s, d) clockwise.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI >
+      IGL_INLINE
+      void order_facets_around_edge(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          size_t s, size_t d, 
+          const std::vector<int>& adj_faces,
+          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
-        // order[0] is the index into adj_face that is immediately after the
-        // pivot face (s, d, pivot point) in clockwise order.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedI>
-        IGL_INLINE
-        void order_facets_around_edge(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                size_t s, size_t d, 
-                const std::vector<int>& adj_faces,
-                const Eigen::PlainObjectBase<DerivedV>& pivot_point,
-                Eigen::PlainObjectBase<DerivedI>& order);
+      // 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
+      // order[0] is the index into adj_face that is immediately after the
+      // pivot face (s, d, pivot point) in clockwise order.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI>
+      IGL_INLINE
+      void order_facets_around_edge(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          size_t s, size_t d, 
+          const std::vector<int>& adj_faces,
+          const Eigen::PlainObjectBase<DerivedV>& pivot_point,
+          Eigen::PlainObjectBase<DerivedI>& order);
     }
   }
 }

include/igl/resolve_duplicated_faces.cpp

@@ -0,0 +1,89 @@
+// 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 "resolve_duplicated_faces.h"
+
+#include "slice.h"
+#include "unique_simplices.h"
+
+template<
+  typename DerivedF1,
+  typename DerivedF2,
+  typename DerivedJ >
+IGL_INLINE void igl::resolve_duplicated_faces(
+    const Eigen::PlainObjectBase<DerivedF1>& F1,
+    Eigen::PlainObjectBase<DerivedF2>& F2,
+    Eigen::PlainObjectBase<DerivedJ>& J) {
+
+  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 on top 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();
+  J.resize(num_kept, 1);
+  std::copy(kept_faces.begin(), kept_faces.end(), J.data());
+  igl::slice(F1, J, 1, F2);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::resolve_duplicated_faces<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::resolve_duplicated_faces<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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/resolve_duplicated_faces.h

@@ -0,0 +1,52 @@
+// 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_COPYLEFT_RESOLVE_DUPLICATED_FACES
+#define IGL_COPYLEFT_RESOLVE_DUPLICATED_FACES
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+
+  // Resolve duplicated faces according to the following rules per unique face:
+  //
+  // 1. If the number of positively oriented faces equals the number of
+  //    negatively oriented faces, remove all duplicated faces at this triangle.
+  // 2. If the number of positively oriented faces equals the number of
+  //    negatively oriented faces plus 1, keeps one of the positively oriented
+  //    face.
+  // 3. If the number of positively oriented faces equals the number of
+  //    negatively oriented faces minus 1, keeps one of the negatively oriented
+  //    face.
+  // 4. If the number of postively oriented faces differ with the number of
+  //    negativley oriented faces by more than 1, the mesh is not orientable.
+  //    An exception will be thrown.
+  //
+  // Inputs:
+  //   F1  #F1 by 3 array of input faces.
+  //
+  // Outputs:
+  //   F2  #F2 by 3 array of output faces without duplicated faces.
+  //   J   #F2 list of indices into F1.
+  template<
+    typename DerivedF1,
+    typename DerivedF2,
+    typename DerivedJ >
+  IGL_INLINE void resolve_duplicated_faces(
+      const Eigen::PlainObjectBase<DerivedF1>& F1,
+      Eigen::PlainObjectBase<DerivedF2>& F2,
+      Eigen::PlainObjectBase<DerivedJ>& J);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "resolve_duplicated_faces.cpp"
+#endif
+
+#endif

include/igl/slice.cpp

@@ -266,4 +266,5 @@ template void igl::slice<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<fl
 template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice<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> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
+template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<long, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

tutorial/CMakeLists.txt

@@ -60,10 +60,6 @@ ENDIF(MSVC)
 
 option(ENABLE_TUTORIALS " " OFF)
 
-add_subdirectory("../external/embree/" "embree")
-include_directories("../external/embree/include")
-list(APPEND SHARED_LIBRARIES "embree")
-
 # endif(True)
 
 # if(True)
@@ -194,8 +190,10 @@ endif(MATLAB_FOUND)
 add_subdirectory("604_Triangle")
 add_subdirectory("605_Tetgen")
 
-add_subdirectory("606_AmbientOcclusion")
-add_subdirectory("607_Picking")
+if (EMBREE_FOUND)
+  add_subdirectory("606_AmbientOcclusion")
+  add_subdirectory("607_Picking")
+endif (EMBREE_FOUND)
 
 if(LIM_FOUND)
   add_subdirectory("608_LIM")