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

Former-commit-id: c0db898b0269321a41c08e6692a75995b2ddd4c8

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -12,8 +12,12 @@
 #include "../cgal/assign_scalar.h"
 #include "../cgal/propagate_winding_numbers.h"
 #include "../cgal/remesh_self_intersections.h"
+#include "../cgal/relabel_small_immersed_cells.h"
+#include "../cgal/extract_cells.h"
+#include "../../extract_manifold_patches.h"
 #include "../../remove_unreferenced.h"
 #include "../../unique_simplices.h"
+#include "../../unique_edge_map.h"
 #include "../../slice.h"
 #include "../../resolve_duplicated_faces.h"
 #include "../../get_seconds.h"
@@ -23,6 +27,7 @@
 
 //#define MESH_BOOLEAN_TIMING
 //#define DOUBLE_CHECK_EXACT_OUTPUT
+//#define SMALL_CELL_REMOVAL
 
 template <
   typename DerivedVA,
@@ -91,6 +96,28 @@ IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
   log_time("resolve_self_intersection");
 #endif
 
+  // Compute edges.
+  Eigen::MatrixXi E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+  // Compute patches
+  Eigen::VectorXi P;
+  const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("patch_extraction");
+#endif
+
+  // Compute cells.
+  Eigen::MatrixXi per_patch_cells;
+  const size_t num_cells =
+    igl::copyleft::cgal::extract_cells(
+        V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("cell_extraction");
+#endif
+
   // Compute winding numbers on each side of each facet.
   const size_t num_faces = F.rows();
   Eigen::MatrixXi W;
@@ -101,7 +128,8 @@ IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
   if (num_faces > 0) 
   {
     valid = valid & 
-      igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
+      igl::copyleft::cgal::propagate_winding_numbers(
+          V, F, uE, uE2E, num_patches, P, num_cells, per_patch_cells, labels, W);
   } else 
   {
     W.resize(0, 4);
@@ -134,6 +162,11 @@ IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
   log_time("compute_output_winding_number");
 #endif
 
+#ifdef SMALL_CELL_REMOVAL
+  igl::copyleft::cgal::relabel_small_immersed_cells(V, F,
+          num_patches, P, num_cells, per_patch_cells, 1e-3, Wr);
+#endif
+
   // Extract boundary separating inside from outside.
   auto index_to_signed_index = [&](size_t i, bool ori) -> int
   {

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -393,10 +393,13 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
-  process_intersecting_boxes();
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("box_intersection_d");
 #endif
+  process_intersecting_boxes();
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  log_time("resolve_intersection");
+#endif
 
   // Convert lIF to Eigen matrix
   assert(lIF.size()%2 == 0);

include/igl/copyleft/cgal/cell_adjacency.cpp

@@ -0,0 +1,31 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "cell_adjacency.h"
+
+template <typename DerivedC>
+IGL_INLINE void igl::copyleft::cgal::cell_adjacency(
+    const Eigen::PlainObjectBase<DerivedC>& per_patch_cells,
+    const size_t num_cells,
+    std::vector<std::set<std::tuple<typename DerivedC::Scalar, bool, size_t> > >&
+    adjacency_list) {
+
+  const size_t num_patches = per_patch_cells.rows();
+  adjacency_list.resize(num_cells);
+  for (size_t i=0; i<num_patches; i++) {
+    const int positive_cell = per_patch_cells(i,0);
+    const int negative_cell = per_patch_cells(i,1);
+    adjacency_list[positive_cell].emplace(negative_cell, false, i);
+    adjacency_list[negative_cell].emplace(positive_cell, true, i);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::copyleft::cgal::cell_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, unsigned long, std::__1::vector<std::__1::set<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::__1::less<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::__1::allocator<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > >, std::__1::allocator<std::__1::set<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::__1::less<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::__1::allocator<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > > > >&);
+#endif

include/igl/copyleft/cgal/cell_adjacency.h

@@ -0,0 +1,50 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_CGAL_CELL_ADJACENCY_H
+#define IGL_COPYLEFT_CGAL_CELL_ADJACENCY_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <set>
+#include <tuple>
+#include <vector>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   per_patch_cells  #P by 2 list of cell labels on each side of each
+      //                    patch.  Cell labels are assumed to be continuous
+      //                    from 0 to #C.
+      //   num_cells        number of cells.
+      //
+      // Outputs:
+      //   adjacency_list  #C array of list of adjcent cell informations.  If
+      //                   cell i and cell j are adjacent via patch x, where i
+      //                   is on the positive side of x, and j is on the
+      //                   negative side.  Then,
+      //                   adjacency_list[i] will contain the entry {j, false, x}
+      //                   and
+      //                   adjacency_list[j] will contain the entry {i, true, x}
+      template < typename DerivedC >
+      IGL_INLINE void cell_adjacency(
+          const Eigen::PlainObjectBase<DerivedC>& per_patch_cells,
+          const size_t num_cells,
+          std::vector<std::set<std::tuple<typename DerivedC::Scalar, bool, size_t> > >&
+          adjacency_list);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cell_adjacency.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/extract_cells.cpp

@@ -154,12 +154,6 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
   std::vector< std::vector<Kernel::Triangle_3 > > 
     triangle_lists(num_components);
   std::vector<std::vector<bool> > in_Is(num_components);
-  for (size_t i=0; i<num_components; i++)
-  {
-    Is[i].resize(components[i].size());
-    std::copy(components[i].begin(), components[i].end(),Is[i].data());
-    submesh_aabb_tree(V,F,Is[i],trees[i],triangle_lists[i],in_Is[i]);
-  }
 
   // Find outer facets, their orientations and cells for each component
   Eigen::VectorXi outer_facets(num_components);
@@ -167,6 +161,8 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
   Eigen::VectorXi outer_cells(num_components);
   for (size_t i=0; i<num_components; i++)
   {
+    Is[i].resize(components[i].size());
+    std::copy(components[i].begin(), components[i].end(),Is[i].data());
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, Is[i], outer_facets[i], flipped);
     outer_facet_orientation[i] = flipped?1:0;
@@ -194,9 +190,9 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     // construct bounding boxes for each component
     DerivedV bbox_min(num_components, 3);
     DerivedV bbox_max(num_components, 3);
-    // Why not just initialize to numeric_limits::min, numeric_limits::max?
-    bbox_min.rowwise() = V.colwise().maxCoeff().eval();
-    bbox_max.rowwise() = V.colwise().minCoeff().eval();
+    // Assuming our mesh (in exact numbers) fits in the range of double.
+    bbox_min.setConstant(std::numeric_limits<double>::max());
+    bbox_max.setConstant(std::numeric_limits<double>::min());
     // Loop over faces
     for (size_t i=0; i<num_faces; i++)
     {
@@ -240,6 +236,9 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
       const size_t num_candidate_comps = candidate_comps.size();
       if (num_candidate_comps == 0) continue;
 
+      // Build aabb tree for this component.
+      submesh_aabb_tree(V,F,Is[i],trees[i],triangle_lists[i],in_Is[i]);
+
       // Get query points on each candidate component: barycenter of
       // outer-facet 
       DerivedV queries(num_candidate_comps, 3);

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -20,6 +20,7 @@
 #include "closest_facet.h"
 #include "assign_scalar.h"
 #include "extract_cells.h"
+#include "cell_adjacency.h"
 
 #include <stdexcept>
 #include <limits>
@@ -54,19 +55,11 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
   };
   tictoc();
 #endif
-  const size_t num_faces = F.rows();
-  //typedef typename DerivedF::Scalar Index;
 
   Eigen::MatrixXi E, uE;
   Eigen::VectorXi EMAP;
   std::vector<std::vector<size_t> > uE2E;
   igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-  bool valid = true;
-  if (!piecewise_constant_winding_number(F, uE, uE2E)) 
-  {
-    std::cerr << "Input mesh is not orientable!" << std::endl;
-    valid = false;
-  }
 
   Eigen::VectorXi P;
   const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
@@ -79,14 +72,61 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
   log_time("cell_extraction");
 #endif
 
-  typedef std::tuple<size_t, bool, size_t> CellConnection;
-  std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
-  for (size_t i=0; i<num_patches; i++) {
-    const int positive_cell = per_patch_cells(i,0);
-    const int negative_cell = per_patch_cells(i,1);
-    cell_adjacency[positive_cell].emplace(negative_cell, false, i);
-    cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+  return propagate_winding_numbers(V, F,
+          uE, uE2E,
+          num_patches, P,
+          num_cells, per_patch_cells,
+          labels, W);
+}
+
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DeriveduE,
+  typename uE2EType,
+  typename DerivedP,
+  typename DerivedC,
+  typename DerivedL,
+  typename DerivedW>
+IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DeriveduE>& uE,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const size_t num_patches,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const size_t num_cells,
+    const Eigen::PlainObjectBase<DerivedC>& C,
+    const Eigen::PlainObjectBase<DerivedL>& labels,
+    Eigen::PlainObjectBase<DerivedW>& W)
+{
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  const auto & tictoc = []() -> double
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  const auto log_time = [&](const std::string& label) -> void {
+    std::cout << "propagate_winding_num." << label << ": "
+      << tictoc() << std::endl;
+  };
+  tictoc();
+#endif
+
+  bool valid = true;
+  if (!piecewise_constant_winding_number(F, uE, uE2E)) 
+  {
+    std::cerr << "Input mesh is not orientable!" << std::endl;
+    valid = false;
   }
+
+  const size_t num_faces = F.rows();
+  typedef std::tuple<typename DerivedC::Scalar, bool, size_t> CellConnection;
+  std::vector<std::set<CellConnection> > cell_adj;
+  igl::copyleft::cgal::cell_adjacency(C, num_cells, cell_adj);
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
   log_time("cell_connectivity");
 #endif
@@ -94,7 +134,7 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
   auto save_cell = [&](const std::string& filename, size_t cell_id) -> void{
     std::vector<size_t> faces;
     for (size_t i=0; i<num_patches; i++) {
-      if ((per_patch_cells.row(i).array() == cell_id).any()) {
+      if ((C.row(i).array() == cell_id).any()) {
         for (size_t j=0; j<num_faces; j++) {
           if ((size_t)P[j] == i) {
             faces.push_back(j);
@@ -140,7 +180,7 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
           size_t curr_idx = Q.front();
           Q.pop();
           int curr_label = cell_labels[curr_idx];
-          for (const auto& neighbor : cell_adjacency[curr_idx]) {
+          for (const auto& neighbor : cell_adj[curr_idx]) {
             if (cell_labels[std::get<0>(neighbor)] == 0) 
             {
               cell_labels[std::get<0>(neighbor)] = curr_label * -1;
@@ -191,7 +231,7 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
 #endif
 
   const size_t outer_patch = P[outer_facet];
-  const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
+  const size_t infinity_cell = C(outer_patch, flipped?1:0);
 
   Eigen::VectorXi patch_labels(num_patches);
   const int INVALID = std::numeric_limits<int>::max();
@@ -214,7 +254,7 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
   while (!Q.empty()) {
     size_t curr_cell = Q.front();
     Q.pop();
-    for (const auto& neighbor : cell_adjacency[curr_cell]) {
+    for (const auto& neighbor : cell_adj[curr_cell]) {
       size_t neighbor_cell, patch_idx;
       bool direction;
       std::tie(neighbor_cell, direction, patch_idx) = neighbor;
@@ -258,8 +298,8 @@ IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
   for (size_t i=0; i<num_faces; i++) 
   {
     const size_t patch = P[i];
-    const size_t positive_cell = per_patch_cells(patch, 0);
-    const size_t negative_cell = per_patch_cells(patch, 1);
+    const size_t positive_cell = C(patch, 0);
+    const size_t negative_cell = C(patch, 1);
     for (size_t j=0; j<num_labels; j++) {
       W(i,j*2  ) = per_cell_W(positive_cell, j);
       W(i,j*2+1) = per_cell_W(negative_cell, j);

include/igl/copyleft/cgal/propagate_winding_numbers.h

@@ -56,6 +56,43 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedL>& labels,
         Eigen::PlainObjectBase<DerivedW>& W);
+
+      // Inputs:
+      //   V  #V by 3 list of vertex positions.
+      //   F  #F by 3 list of triangle indices into V.
+      //   uE    #uE by 2 list of vertex_indices, represents undirected edges.
+      //   uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+      //   num_patches  number of patches
+      //   P  #F list of patch ids.
+      //   num_cells    number of cells
+      //   C  #P by 2 list of cell ids on each side of each patch.
+      //   labels  #F list of facet labels ranging from 0 to k-1.
+      // Output:
+      //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
+      //      number on the positive side of facet ``i`` with respect to the
+      //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
+      //      number on the negative side of facet ``i`` with respect to the
+      //      facets labeled ``j``.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DeriveduE,
+        typename uE2EType,
+        typename DerivedP,
+        typename DerivedC,
+        typename DerivedL,
+        typename DerivedW>
+      IGL_INLINE bool propagate_winding_numbers(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const size_t num_patches,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const size_t num_cells,
+        const Eigen::PlainObjectBase<DerivedC>& C,
+        const Eigen::PlainObjectBase<DerivedL>& labels,
+        Eigen::PlainObjectBase<DerivedW>& W);
     }
   }
 }

include/igl/copyleft/cgal/relabel_small_immersed_cells.cpp

@@ -0,0 +1,121 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "relabel_small_immersed_cells.h"
+#include "../../centroid.h"
+#include "assign_scalar.h"
+#include "cell_adjacency.h"
+
+#include <vector>
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedC,
+  typename FT,
+  typename DerivedW>
+IGL_INLINE void igl::copyleft::cgal::relabel_small_immersed_cells(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const size_t num_patches,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const size_t num_cells,
+    const Eigen::PlainObjectBase<DerivedC>& C,
+    const FT vol_threashold,
+    Eigen::PlainObjectBase<DerivedW>& W)
+{
+  const size_t num_vertices = V.rows();
+  const size_t num_faces = F.rows();
+  typedef std::tuple<typename DerivedC::Scalar, bool, size_t> CellConnection;
+  std::vector<std::set<CellConnection> > cell_adj;
+  igl::copyleft::cgal::cell_adjacency(C, num_cells, cell_adj);
+
+  Eigen::MatrixXd VV(V.rows(), V.cols());
+  for (size_t i=0; i<num_vertices; i++) {
+    igl::copyleft::cgal::assign_scalar(V(i,0), VV(i,0));
+    igl::copyleft::cgal::assign_scalar(V(i,1), VV(i,1));
+    igl::copyleft::cgal::assign_scalar(V(i,2), VV(i,2));
+  }
+
+  auto compute_cell_volume = [&](size_t cell_id) {
+    std::vector<short> is_involved(num_patches, 0);
+    for (size_t i=0; i<num_patches; i++) {
+      if (C(i,0) == cell_id) {
+        // cell is on positive side of patch i.
+        is_involved[i] = 1;
+      }
+      if (C(i,1) == cell_id) {
+        // cell is on negative side of patch i.
+        is_involved[i] = -1;
+      }
+    }
+
+    std::vector<size_t> involved_positive_faces;
+    std::vector<size_t> involved_negative_faces;
+    for (size_t i=0; i<num_faces; i++) {
+      switch (is_involved[P[i]]) {
+        case 1:
+          involved_negative_faces.push_back(i);
+          break;
+        case -1:
+          involved_positive_faces.push_back(i);
+          break;
+      }
+    }
+
+    const size_t num_positive_faces = involved_positive_faces.size();
+    const size_t num_negative_faces = involved_negative_faces.size();
+    DerivedF selected_faces(num_positive_faces + num_negative_faces, 3);
+    for (size_t i=0; i<num_positive_faces; i++) {
+      selected_faces.row(i) = F.row(involved_positive_faces[i]);
+    }
+    for (size_t i=0; i<num_negative_faces; i++) {
+      selected_faces.row(num_positive_faces+i) =
+        F.row(involved_negative_faces[i]).reverse();
+    }
+
+    Eigen::VectorXd c(3);
+    double vol;
+    igl::centroid(VV, selected_faces, c, vol);
+    return vol;
+  };
+
+  std::vector<typename DerivedV::Scalar> cell_volumes(num_cells);
+  for (size_t i=0; i<num_cells; i++) {
+    cell_volumes[i] = compute_cell_volume(i);
+  }
+
+  std::vector<typename DerivedW::Scalar> cell_values(num_cells);
+  for (size_t i=0; i<num_faces; i++) {
+    cell_values[C(P[i], 0)] = W(i, 0);
+    cell_values[C(P[i], 1)] = W(i, 1);
+  }
+
+  for (size_t i=1; i<num_cells; i++) {
+    std::cout << cell_volumes[i] << std::endl;
+    if (cell_volumes[i] >= vol_threashold) continue;
+    std::set<typename DerivedW::Scalar> neighbor_values;
+    const auto neighbors = cell_adj[i];
+    for (const auto& entry : neighbors) {
+      const auto& j = std::get<0>(entry);
+      neighbor_values.insert(cell_values[j]);
+    }
+    // If cell i is immersed, assign its value to be the immersed value.
+    if (neighbor_values.size() == 1) {
+      cell_values[i] = *neighbor_values.begin();
+    }
+  }
+
+  for (size_t i=0; i<num_faces; i++) {
+    W(i,0) = cell_values[C(P[i], 0)];
+    W(i,1) = cell_values[C(P[i], 1)];
+  }
+}
+

include/igl/copyleft/cgal/relabel_small_immersed_cells.h

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_RELABEL_SMALL_IMMERSED_CELLS
+#define IGL_RELABEL_SMALL_IMMERSED_CELLS
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   V  #V by 3 list of vertex positions.
+      //   F  #F by 3 list of triangle indices into V.
+      //   num_patches  number of patches
+      //   P  #F list of patch ids.
+      //   num_cells    number of cells
+      //   C  #P by 2 list of cell ids on each side of each patch.
+      //   vol_threshold  Volume threshold, cells smaller than this
+      //                  and is completely immersed will be relabeled.
+      //
+      // In/Output:
+      //   W  #F by 2 cell labels.  W(i,0) is the label on the positive side of
+      //      face i, W(i,1) is the label on the negative side of face i.  W
+      //      will be modified in place by this method.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename DerivedC,
+        typename FT,
+        typename DerivedW>
+      IGL_INLINE void relabel_small_immersed_cells(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const size_t num_patches,
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          const size_t num_cells,
+          const Eigen::PlainObjectBase<DerivedC>& C,
+          const FT vol_threashold,
+          Eigen::PlainObjectBase<DerivedW>& W);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "relabel_small_immersed_cells.cpp"
+#endif
+#endif