Replace old mesh boolean with the beta one.

Former-commit-id: c0d8135b9ae6e7f682d248dc1f9563b5cdefd6fb

include/igl/boolean/mesh_boolean.cpp

@@ -1,379 +1,413 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
 // 
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
+//
 #include "mesh_boolean.h"
 #include <igl/cgal/assign_scalar.h>
-#include <igl/per_face_normals.h>
-#include <igl/boundary_facets.h>
-#include <igl/exterior_edges.h>
-#include <igl/cgal/peel_outer_hull_layers.h>
+#include <igl/cgal/propagate_winding_numbers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
-#include <igl/mod.h>
-#include <igl/unique_simplices.h>
+
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <iostream>
 
-//#define IGL_MESH_BOOLEAN_DEBUG
+#include "../writePLY.h"
+#include "../writeDMAT.h"
+
+namespace igl {
+    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::cgal::RemeshSelfIntersectionsParam params;
+
+                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());
+
+                // 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;
+                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]);
+                }
+            }
+
+            typedef Eigen::Matrix<int, 1, Eigen::Dynamic> WindingNumbers;
+            typedef std::function<int(const WindingNumbers&)> WindingNumberOperation;
+
+            WindingNumberOperation binary_union() {
+                return [](const WindingNumbers& win_nums) -> int{
+                    return win_nums[0] > 0 || win_nums[1] > 0;
+                };
+            }
+
+            WindingNumberOperation binary_intersect() {
+                return [](const WindingNumbers& win_nums) -> int{
+                    return win_nums[0] > 0 && win_nums[1] > 0;
+                };
+            }
+
+            WindingNumberOperation binary_difference() {
+                return [](const WindingNumbers& win_nums) -> int{
+                    return win_nums[0] > 0 && win_nums[1] <= 0;
+                };
+            }
+
+            WindingNumberOperation binary_xor() {
+                return [](const WindingNumbers& win_nums) -> int{
+                    return (win_nums[0] > 0 && win_nums[1] <= 0) ||
+                           (win_nums[0] <= 0 && win_nums[1] > 0);
+                };
+            }
+
+            WindingNumberOperation binary_resolve() {
+                return [](const WindingNumbers& win_nums) -> int{
+                    return true;
+                };
+            }
+
+            typedef std::function<short(int, int)> ToKeepFunc;
+            ToKeepFunc keep_inside() {
+                return [](int out_w, int in_w) -> short {
+                    if (in_w > 0 && out_w <= 0) return 1;
+                    else if (in_w <= 0 && out_w > 0) return -1;
+                    else return 0;
+                };
+            }
+
+            ToKeepFunc keep_all() {
+                return [](int out_w, int in_w) -> short {
+                    return true;
+                };
+            }
+        }
+    }
+}
+
 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);});
+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_beta(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);
     }
-#ifdef IGL_MESH_BOOLEAN_DEBUG
-    cout<<"#F:  "<<F.rows()<<endl;
-    cout<<"#CF: "<<CF.rows()<<endl;
-#endif
-  };
 
-#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));
-      }
+    // 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);
     }
-  }
-
-  if(type == MESH_BOOLEAN_TYPE_RESOLVE)
-  {
-    FC = CF;
-    VC = CV;
-    J = CJ;
-    I = CI;
-    return;
-  }
-
-#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;
-  };
-
-#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);
+
+    // 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));
         }
-        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(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)++;
-    }
-    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 ?");
-      }
     }
-    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]);
+
+    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;
+    WindingNumberOperation op;
+    ToKeepFunc keep;
+    switch (type) {
+        case MESH_BOOLEAN_TYPE_UNION:
+            op = binary_union();
+            keep = keep_inside();
+            break;
+        case MESH_BOOLEAN_TYPE_INTERSECT:
+            op = binary_intersect();
+            keep = keep_inside();
+            break;
+        case MESH_BOOLEAN_TYPE_MINUS:
+            op = binary_difference();
+            keep = keep_inside();
+            break;
+        case MESH_BOOLEAN_TYPE_XOR:
+            op = binary_xor();
+            keep = keep_inside();
+            break;
+        case MESH_BOOLEAN_TYPE_RESOLVE:
+            op = binary_resolve();
+            keep = keep_all();
+            break;
+        default:
+            throw std::runtime_error("Unsupported boolean type.");
+    }
+
+    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>;
+
+    igl::boolean::per_face_winding_number_binary_operation(
+            VA, FA, VB, FB, op, keep, resolve_func, VC, FC, J);
 }
 
 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
@@ -381,24 +415,4 @@ IGL_INLINE void igl::boolean::mesh_boolean(
 template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::boolean::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-// This is a hack to discuss. I'm not sure why this _doesn't_ create
-// duplicate symbols.
-#include <igl/remove_unreferenced.cpp>
-template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-#include <igl/cgal/peel_outer_hull_layers.cpp>
-template unsigned long
-igl::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-#include <igl/cgal/outer_hull.cpp>
-template void igl::cgal::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-#include <igl/slice.cpp>
-#include <igl/barycenter.cpp>
-template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
-#include <igl/mod.cpp>
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > igl::mod<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int);
-#include <igl/outer_element.cpp>
-template void igl::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-#include <igl/colon.cpp>
-template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
-
-
 #endif

include/igl/boolean/mesh_boolean.h

@@ -1,10 +1,11 @@
 // 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_BOOLEAN_MESH_BOOLEAN_H
 #define IGL_BOOLEAN_MESH_BOOLEAN_H
 
@@ -17,54 +18,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 +48,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 +57,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/boolean/mesh_boolean_beta.cpp

@@ -1,411 +0,0 @@
-// 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 "mesh_boolean_beta.h"
-#include <igl/cgal/assign_scalar.h>
-#include <igl/cgal/propagate_winding_numbers.h>
-#include <igl/cgal/remesh_self_intersections.h>
-#include <igl/remove_unreferenced.h>
-
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-
-#include "../writePLY.h"
-#include "../writeDMAT.h"
-
-namespace igl {
-    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::cgal::RemeshSelfIntersectionsParam params;
-
-                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());
-
-                // 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;
-                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]);
-                }
-            }
-
-            typedef Eigen::Matrix<int, 1, Eigen::Dynamic> WindingNumbers;
-            typedef std::function<int(const WindingNumbers&)> WindingNumberOperation;
-
-            WindingNumberOperation binary_union() {
-                return [](const WindingNumbers& win_nums) -> int{
-                    return win_nums[0] > 0 || win_nums[1] > 0;
-                };
-            }
-
-            WindingNumberOperation binary_intersect() {
-                return [](const WindingNumbers& win_nums) -> int{
-                    return win_nums[0] > 0 && win_nums[1] > 0;
-                };
-            }
-
-            WindingNumberOperation binary_difference() {
-                return [](const WindingNumbers& win_nums) -> int{
-                    return win_nums[0] > 0 && win_nums[1] <= 0;
-                };
-            }
-
-            WindingNumberOperation binary_xor() {
-                return [](const WindingNumbers& win_nums) -> int{
-                    return (win_nums[0] > 0 && win_nums[1] <= 0) ||
-                           (win_nums[0] <= 0 && win_nums[1] > 0);
-                };
-            }
-
-            WindingNumberOperation binary_resolve() {
-                return [](const WindingNumbers& win_nums) -> int{
-                    return true;
-                };
-            }
-
-            typedef std::function<short(int, int)> ToKeepFunc;
-            ToKeepFunc keep_inside() {
-                return [](int out_w, int in_w) -> short {
-                    if (in_w > 0 && out_w <= 0) return 1;
-                    else if (in_w <= 0 && out_w > 0) return -1;
-                    else return 0;
-                };
-            }
-
-            ToKeepFunc keep_all() {
-                return [](int out_w, int in_w) -> short {
-                    return true;
-                };
-            }
-        }
-    }
-}
-
-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_beta(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);
-    }
-
-    // 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);
-    }
-
-    // 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];
-    }
-
-
-    // Finally, remove duplicated faces and unreferenced vertices.
-    {
-        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);
-    }
-}
-
-template <
-typename DerivedVA,
-typename DerivedFA,
-typename DerivedVB,
-typename DerivedFB,
-typename DerivedVC,
-typename DerivedFC,
-typename DerivedJ>
-IGL_INLINE void igl::boolean::mesh_boolean_beta(
-        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;
-    WindingNumberOperation op;
-    ToKeepFunc keep;
-    switch (type) {
-        case MESH_BOOLEAN_TYPE_UNION:
-            op = binary_union();
-            keep = keep_inside();
-            break;
-        case MESH_BOOLEAN_TYPE_INTERSECT:
-            op = binary_intersect();
-            keep = keep_inside();
-            break;
-        case MESH_BOOLEAN_TYPE_MINUS:
-            op = binary_difference();
-            keep = keep_inside();
-            break;
-        case MESH_BOOLEAN_TYPE_XOR:
-            op = binary_xor();
-            keep = keep_inside();
-            break;
-        case MESH_BOOLEAN_TYPE_RESOLVE:
-            op = binary_resolve();
-            keep = keep_all();
-            break;
-        default:
-            throw std::runtime_error("Unsupported boolean type.");
-    }
-
-    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>;
-
-    igl::boolean::per_face_winding_number_binary_operation(
-            VA, FA, VB, FB, op, keep, resolve_func, VC, FC, J);
-}
-
-template <
-typename DerivedVA,
-typename DerivedFA,
-typename DerivedVB,
-typename DerivedFB,
-typename DerivedVC,
-typename DerivedFC>
-IGL_INLINE void igl::boolean::mesh_boolean_beta(
-        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_beta(VA,FA,VB,FB,type,VC,FC,J);
-}

include/igl/boolean/mesh_boolean_beta.h

@@ -1,85 +0,0 @@
-// 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_BOOLEAN_MESH_BOOLEAN_BETA_H
-#define IGL_BOOLEAN_MESH_BOOLEAN_BETA_H
-
-#include <igl/igl_inline.h>
-#include "MeshBooleanType.h"
-#include <Eigen/Core>
-#include <functional>
-
-namespace igl
-{
-  namespace boolean
-  {
-      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,
-      typename DerivedVB,
-      typename DerivedFB,
-      typename DerivedVC,
-      typename DerivedFC,
-      typename DerivedJ>
-    IGL_INLINE void mesh_boolean_beta(
-      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);
-
-    template <
-      typename DerivedVA,
-      typename DerivedFA,
-      typename DerivedVB,
-      typename DerivedFB,
-      typename DerivedVC,
-      typename DerivedFC>
-    IGL_INLINE void mesh_boolean_beta(
-      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);
-
-  }
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "mesh_boolean_beta.cpp"
-#endif
-
-#endif