Merge branch 'master' of github.com:libigl/libigl

Former-commit-id: 9c03252602e589d47f07cb22ada79d2ea6f01954

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -354,7 +354,10 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     tit != T.end(); 
     ++tit)
   {
-    boxes.push_back(Box(tit->bbox(), tit));
+    if (!tit->is_degenerate())
+    {
+      boxes.push_back(Box(tit->bbox(), tit));
+    }
   }
   // Leapfrog callback
   std::function<void(const Box &a,const Box &b)> cb = 
@@ -831,10 +834,10 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   Index fb = b.handle()-T.begin();
   const Triangle_3 & A = *a.handle();
   const Triangle_3 & B = *b.handle();
-  // I'm not going to deal with degenerate triangles, though at some point we
-  // should
-  assert(!a.handle()->is_degenerate());
-  assert(!b.handle()->is_degenerate());
+  //// I'm not going to deal with degenerate triangles, though at some point we
+  //// should
+  //assert(!a.handle()->is_degenerate());
+  //assert(!b.handle()->is_degenerate());
   // Number of combinatorially shared vertices
   Index comb_shared_vertices = 0;
   // Number of geometrically shared vertices (*not* including combinatorially

include/igl/copyleft/cgal/closest_facet.cpp

@@ -10,6 +10,7 @@
 
 #include <vector>
 #include <stdexcept>
+#include <unordered_map>
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
@@ -18,343 +19,412 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 #include "order_facets_around_edge.h"
+#include "../../vertex_triangle_adjacency.h"
+#include "../../writePLY.h"
 
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedI,
-    typename DerivedP,
-    typename DerivedR,
-    typename DerivedS >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedR,
+  typename DerivedS >
 IGL_INLINE void igl::copyleft::cgal::closest_facet(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedI>& I,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        Eigen::PlainObjectBase<DerivedR>& R,
-        Eigen::PlainObjectBase<DerivedS>& S) {
-    typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    typedef Kernel::Point_3 Point_3;
-    typedef Kernel::Plane_3 Plane_3;
-    typedef Kernel::Segment_3 Segment_3;
-    typedef Kernel::Triangle_3 Triangle;
-    typedef std::vector<Triangle>::iterator Iterator;
-    typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-    typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-    typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-
-    if (F.rows() <= 0 || I.rows() <= 0) {
-        throw std::runtime_error(
-                "Closest facet cannot be computed on empty mesh.");
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedI>& I,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S) {
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+  typedef Kernel::Point_3 Point_3;
+  typedef Kernel::Plane_3 Plane_3;
+  typedef Kernel::Segment_3 Segment_3;
+  typedef Kernel::Triangle_3 Triangle;
+  typedef std::vector<Triangle>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+  if (F.rows() <= 0 || I.rows() <= 0) {
+    throw std::runtime_error(
+        "Closest facet cannot be computed on empty mesh.");
+  }
+
+  std::vector<std::vector<size_t> > VF;
+  std::vector<std::vector<size_t> > VFi;
+  igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
+
+  std::vector<bool> in_I(F.rows(), false);
+  const size_t num_faces = I.rows();
+  std::vector<Triangle> triangles;
+  for (size_t i=0; i<num_faces; i++) {
+    const Eigen::Vector3i f = F.row(I(i, 0));
+    in_I[I(i,0)] = true;
+    triangles.emplace_back(
+        Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+        Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+        Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+    if (triangles.back().is_degenerate()) {
+      throw std::runtime_error(
+          "Input facet components contains degenerated triangles");
+    }
+  }
+  Tree tree(triangles.begin(), triangles.end());
+  tree.accelerate_distance_queries();
+
+  auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
+    const auto& f = F.row(fid).eval();
+    Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+    Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+    Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+    auto ori = CGAL::orientation(v0, v1, v2, p);
+    switch (ori) {
+      case CGAL::POSITIVE:
+        return true;
+      case CGAL::NEGATIVE:
+        return false;
+      case CGAL::COPLANAR:
+        // Warning:
+        // This can only happen if fid contains a boundary edge.
+        // Catergorized this ambiguous case as negative side.
+        return false;
+      default:
+        throw std::runtime_error("Unknown CGAL state.");
+    }
+    return false;
+  };
+
+  auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
+    const auto& f = F.row(fid);
+    if      ((size_t)f[0] == s && (size_t)f[1] == d) return false;
+    else if ((size_t)f[1] == s && (size_t)f[2] == d) return false;
+    else if ((size_t)f[2] == s && (size_t)f[0] == d) return false;
+    else if ((size_t)f[0] == d && (size_t)f[1] == s) return true;
+    else if ((size_t)f[1] == d && (size_t)f[2] == s) return true;
+    else if ((size_t)f[2] == d && (size_t)f[0] == s) return true;
+    else {
+      throw std::runtime_error(
+          "Cannot compute orientation due to incorrect connectivity");
+      return false;
+    }
+  };
+  auto index_to_signed_index = [&](size_t index, bool ori) -> int{
+    return (index+1) * (ori? 1:-1);
+  };
+  //auto signed_index_to_index = [&](int signed_index) -> size_t {
+  //    return abs(signed_index) - 1;
+  //};
+
+  enum ElementType { VERTEX, EDGE, FACE };
+  auto determine_element_type = [&](const Point_3& p, const size_t fid,
+      size_t& element_index) {
+    const auto& tri = triangles[fid];
+    const Point_3 p0 = tri[0];
+    const Point_3 p1 = tri[1];
+    const Point_3 p2 = tri[2];
+
+    if (p == p0) { element_index = 0; return VERTEX; }
+    if (p == p1) { element_index = 1; return VERTEX; }
+    if (p == p2) { element_index = 2; return VERTEX; }
+    if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+    if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+    if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+    element_index = 0;
+    return FACE;
+  };
+
+  auto process_edge_case = [&](
+      size_t query_idx,
+      const size_t s, const size_t d,
+      size_t preferred_facet,
+      bool& orientation) {
+    Point_3 query_point(
+        P(query_idx, 0),
+        P(query_idx, 1),
+        P(query_idx, 2));
+
+    size_t corner_idx = std::numeric_limits<size_t>::max();
+    if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 1)) ||
+        (s == F(preferred_facet, 1) && d == F(preferred_facet, 0))) {
+      corner_idx = 2;
+    } else if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 2)) ||
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 0))) {
+      corner_idx = 1;
+    } else if ((s == F(preferred_facet, 1) && d == F(preferred_facet, 2)) ||
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 1))) {
+      corner_idx = 0;
+    } else {
+      std::cerr << "s: " << s << "\t d:" << d << std::endl;
+      std::cerr << F.row(preferred_facet) << std::endl;
+      throw std::runtime_error(
+          "Invalid connectivity, edge does not belong to facet");
     }
 
-    const size_t num_faces = I.rows();
-    std::vector<Triangle> triangles;
-    for (size_t i=0; i<num_faces; i++) {
-        const Eigen::Vector3i f = F.row(I(i, 0));
-        triangles.emplace_back(
-                Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
-                Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
-                Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
-        if (triangles.back().is_degenerate()) {
-            throw std::runtime_error(
-                    "Input facet components contains degenerated triangles");
-        }
+    auto ueid = EMAP(preferred_facet + corner_idx * F.rows());
+    auto eids = uE2E[ueid];
+    std::vector<size_t> intersected_face_indices;
+    for (auto eid : eids) {
+      const size_t fid = eid % F.rows();
+      if (in_I[fid]) {
+        intersected_face_indices.push_back(fid);
+      }
     }
-    Tree tree(triangles.begin(), triangles.end());
-    tree.accelerate_distance_queries();
-
-    auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
-        const auto& f = F.row(fid).eval();
-        Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
-        Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
-        Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
-        auto ori = CGAL::orientation(v0, v1, v2, p);
-        switch (ori) {
-            case CGAL::POSITIVE:
-                return true;
-            case CGAL::NEGATIVE:
-                return false;
-            case CGAL::COPLANAR:
-                throw std::runtime_error(
-                        "It seems input mesh contains self intersection");
-            default:
-                throw std::runtime_error("Unknown CGAL state.");
-        }
-        return false;
-    };
 
-    auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
-        const auto& f = F.row(fid);
-        if      ((size_t)f[0] == s && (size_t)f[1] == d) return false;
-        else if ((size_t)f[1] == s && (size_t)f[2] == d) return false;
-        else if ((size_t)f[2] == s && (size_t)f[0] == d) return false;
-        else if ((size_t)f[0] == d && (size_t)f[1] == s) return true;
-        else if ((size_t)f[1] == d && (size_t)f[2] == s) return true;
-        else if ((size_t)f[2] == d && (size_t)f[0] == s) return true;
-        else {
-            throw std::runtime_error(
-                    "Cannot compute orientation due to incorrect connectivity");
-            return false;
-        }
-    };
-    auto index_to_signed_index = [&](size_t index, bool ori) -> int{
-        return (index+1) * (ori? 1:-1);
-    };
-    //auto signed_index_to_index = [&](int signed_index) -> size_t {
-    //    return abs(signed_index) - 1;
-    //};
-
-    enum ElementType { VERTEX, EDGE, FACE };
-    auto determine_element_type = [&](const Point_3& p, const size_t fid,
-            size_t& element_index) {
-        const auto& tri = triangles[fid];
-        const Point_3 p0 = tri[0];
-        const Point_3 p1 = tri[1];
-        const Point_3 p2 = tri[2];
-
-        if (p == p0) { element_index = 0; return VERTEX; }
-        if (p == p1) { element_index = 1; return VERTEX; }
-        if (p == p2) { element_index = 2; return VERTEX; }
-        if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
-        if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
-        if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
-
-        element_index = 0;
-        return FACE;
-    };
+    const size_t num_intersected_faces = intersected_face_indices.size();
+    std::vector<int> intersected_face_signed_indices(num_intersected_faces);
+    std::transform(
+        intersected_face_indices.begin(),
+        intersected_face_indices.end(),
+        intersected_face_signed_indices.begin(),
+        [&](size_t index) {
+        return index_to_signed_index(
+          index, get_orientation(index, s,d));
+        });
+
+    assert(num_intersected_faces >= 1);
+    if (num_intersected_faces == 1) {
+      // The edge must be a boundary edge.  Thus, the orientation can be
+      // simply determined by checking if the query point is on the
+      // positive side of the facet.
+      const size_t fid = intersected_face_indices[0];
+      orientation = on_the_positive_side(fid, query_point);
+      return fid;
+    }
 
-    auto process_edge_case = [&](
-            size_t query_idx,
-            const size_t s, const size_t d,
-            size_t preferred_facet,
-            bool& orientation) {
-
-        Point_3 mid_edge_point(
-                (V(s,0) + V(d,0)) * 0.5,
-                (V(s,1) + V(d,1)) * 0.5,
-                (V(s,2) + V(d,2)) * 0.5);
-        Point_3 query_point(
-                P(query_idx, 0),
-                P(query_idx, 1),
-                P(query_idx, 2));
-
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(mid_edge_point, query_point),
-                std::back_inserter(intersected_faces));
-
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::vector<int> intersected_face_signed_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-        std::transform(
-                intersected_face_indices.begin(),
-                intersected_face_indices.end(),
-                intersected_face_signed_indices.begin(),
-                [&](size_t index) {
-                    return index_to_signed_index(
-                        index, get_orientation(index, s,d));
-                });
-
-        assert(num_intersected_faces >= 1);
-        if (num_intersected_faces == 1) {
-            // The edge must be a boundary edge.  Thus, the orientation can be
-            // simply determined by checking if the query point is on the
-            // positive side of the facet.
-            const size_t fid = intersected_face_indices[0];
-            orientation = on_the_positive_side(fid, query_point);
-            return fid;
+    Eigen::VectorXi order;
+    DerivedP pivot = P.row(query_idx).eval();
+    igl::copyleft::cgal::order_facets_around_edge(V, F, s, d,
+        intersected_face_signed_indices,
+        pivot, order);
+
+    // Although first and last are equivalent, make the choice based on
+    // preferred_facet.
+    const size_t first = order[0];
+    const size_t last = order[num_intersected_faces-1];
+    if (intersected_face_indices[first] == preferred_facet) {
+      orientation = intersected_face_signed_indices[first] < 0;
+      return intersected_face_indices[first];
+    } else if (intersected_face_indices[last] == preferred_facet) {
+      orientation = intersected_face_signed_indices[last] > 0;
+      return intersected_face_indices[last];
+    } else {
+      orientation = intersected_face_signed_indices[order[0]] < 0;
+      return intersected_face_indices[order[0]];
+    }
+  };
+
+  auto process_face_case = [&](
+      const size_t query_idx, const Point_3& closest_point,
+      const size_t fid, bool& orientation) {
+    const auto& f = F.row(I(fid, 0));
+    return process_edge_case(query_idx, f[0], f[1], I(fid, 0), orientation);
+  };
+
+  auto process_vertex_case = [&](const size_t query_idx, size_t s,
+      size_t preferred_facet, bool& orientation) {
+    const Point_3 query_point(
+        P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
+    const Point_3 closest_point(V(s, 0), V(s, 1), V(s, 2));
+    std::vector<size_t> adj_faces;
+    std::vector<size_t> adj_face_corners;
+    {
+      // Gether adj faces to s within I.
+      const auto& all_adj_faces = VF[s];
+      const auto& all_adj_face_corners = VFi[s];
+      const size_t num_all_adj_faces = all_adj_faces.size();
+      for (size_t i=0; i<num_all_adj_faces; i++) {
+        const size_t fid = all_adj_faces[i];
+        if (in_I[fid]) {
+          adj_faces.push_back(fid);
+          adj_face_corners.push_back(all_adj_face_corners[i]);
         }
+      }
+    }
+    const size_t num_adj_faces = adj_faces.size();
+    assert(num_adj_faces > 0);
+
+    std::set<size_t> adj_vertices_set;
+    std::unordered_multimap<size_t, size_t> v2f;
+    for (size_t i=0; i<num_adj_faces; i++) {
+      const size_t fid = adj_faces[i];
+      const size_t cid = adj_face_corners[i];
+      const auto& f = F.row(adj_faces[i]);
+      const size_t next = f[(cid+1)%3];
+      const size_t prev = f[(cid+2)%3];
+      adj_vertices_set.insert(next);
+      adj_vertices_set.insert(prev);
+      v2f.insert({{next, fid}, {prev, fid}});
+    }
+    const size_t num_adj_vertices = adj_vertices_set.size();
+    std::vector<size_t> adj_vertices(num_adj_vertices);
+    std::copy(adj_vertices_set.begin(), adj_vertices_set.end(),
+        adj_vertices.begin());
+
+    std::vector<Point_3> adj_points;
+    for (size_t vid : adj_vertices) {
+      adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
+    }
 
-        Eigen::VectorXi order;
-        DerivedP pivot = P.row(query_idx).eval();
-        igl::copyleft::cgal::order_facets_around_edge(V, F, s, d,
-                intersected_face_signed_indices,
-                pivot, order);
-
-        // Although first and last are equivalent, make the choice based on
-        // preferred_facet.
-        const size_t first = order[0];
-        const size_t last = order[num_intersected_faces-1];
-        if (intersected_face_indices[first] == preferred_facet) {
-            orientation = intersected_face_signed_indices[first] < 0;
-            return intersected_face_indices[first];
-        } else if (intersected_face_indices[last] == preferred_facet) {
-            orientation = intersected_face_signed_indices[last] > 0;
-            return intersected_face_indices[last];
-        } else {
-            orientation = intersected_face_signed_indices[order[0]] < 0;
-            return intersected_face_indices[order[0]];
+    // A plane is on the exterior if all adj_points lies on or to
+    // one side of the plane.
+    auto is_on_exterior = [&](const Plane_3& separator) {
+      size_t positive=0;
+      size_t negative=0;
+      size_t coplanar=0;
+      for (const auto& point : adj_points) {
+        switch(separator.oriented_side(point)) {
+          case CGAL::ON_POSITIVE_SIDE:
+            positive++;
+            break;
+          case CGAL::ON_NEGATIVE_SIDE:
+            negative++;
+            break;
+          case CGAL::ON_ORIENTED_BOUNDARY:
+            coplanar++;
+            break;
+          default:
+            throw "Unknown plane-point orientation";
         }
+      }
+      auto query_orientation = separator.oriented_side(query_point);
+      if (query_orientation == CGAL::ON_ORIENTED_BOUNDARY &&
+          (positive == 0 && negative == 0)) {
+        // All adj vertices and query point are coplanar.
+        // In this case, all separators are equally valid.
+        return true;
+      } else {
+        bool r = (positive == 0 && query_orientation == CGAL::POSITIVE)
+          || (negative == 0 && query_orientation == CGAL::NEGATIVE);
+        return r;
+      }
     };
 
-    auto process_face_case = [&](
-            const size_t query_idx, const size_t fid, bool& orientation) {
-        const auto& f = F.row(I(fid, 0));
-        return process_edge_case(query_idx, f[0], f[1], I(fid, 0), orientation);
-    };
-
-    auto process_vertex_case = [&](const size_t query_idx, size_t s,
-            size_t preferred_facet, bool& orientation) {
-        Point_3 closest_point(V(s, 0), V(s, 1), V(s, 2));
-        Point_3 query_point(P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
-
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(closest_point, query_point),
-                std::back_inserter(intersected_faces));
-
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-
-        std::set<size_t> adj_vertices_set;
-        for (auto fid : intersected_face_indices) {
-            const auto& f = F.row(fid);
-            if ((size_t)f[0] != s) adj_vertices_set.insert(f[0]);
-            if ((size_t)f[1] != s) adj_vertices_set.insert(f[1]);
-            if ((size_t)f[2] != s) adj_vertices_set.insert(f[2]);
+    size_t d = std::numeric_limits<size_t>::max();
+    for (size_t i=0; i<num_adj_vertices; i++) {
+      const size_t vi = adj_vertices[i];
+      for (size_t j=i+1; j<num_adj_vertices; j++) {
+        Plane_3 separator(closest_point, adj_points[i], adj_points[j]);
+        if (separator.is_degenerate()) {
+          continue;
         }
-        const size_t num_adj_vertices = adj_vertices_set.size();
-        std::vector<size_t> adj_vertices(num_adj_vertices);
-        std::copy(adj_vertices_set.begin(), adj_vertices_set.end(),
-                adj_vertices.begin());
-
-        std::vector<Point_3> adj_points;
-        for (size_t vid : adj_vertices) {
-            adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
+        if (is_on_exterior(separator)) {
+          if (!CGAL::collinear(
+                query_point, adj_points[i], closest_point)) {
+            d = vi;
+            break;
+          } else {
+            d = adj_vertices[j];
+            assert(!CGAL::collinear(
+                  query_point, adj_points[j], closest_point));
+            break;
+          }
         }
-
-        // A plane is on the exterior if all adj_points lies on or to
-        // one side of the plane.
-        auto is_on_exterior = [&](const Plane_3& separator) {
-            size_t positive=0;
-            size_t negative=0;
-            size_t coplanar=0;
-            for (const auto& point : adj_points) {
-                switch(separator.oriented_side(point)) {
-                    case CGAL::ON_POSITIVE_SIDE:
-                        positive++;
-                        break;
-                    case CGAL::ON_NEGATIVE_SIDE:
-                        negative++;
-                        break;
-                    case CGAL::ON_ORIENTED_BOUNDARY:
-                        coplanar++;
-                        break;
-                    default:
-                        throw "Unknown plane-point orientation";
-                }
-            }
-            auto query_orientation = separator.oriented_side(query_point);
-            bool r = (positive == 0 && query_orientation == CGAL::POSITIVE)
-                || (negative == 0 && query_orientation == CGAL::NEGATIVE);
-            return r;
-        };
-
-        size_t d = std::numeric_limits<size_t>::max();
-        for (size_t i=0; i<num_adj_vertices; i++) {
-            const size_t vi = adj_vertices[i];
-            for (size_t j=i+1; j<num_adj_vertices; j++) {
-                Plane_3 separator(closest_point, adj_points[i], adj_points[j]);
-                if (separator.is_degenerate()) {
-                    throw std::runtime_error(
-                            "Input mesh contains degenerated faces");
-                }
-                if (is_on_exterior(separator)) {
-                    d = vi;
-                    assert(!CGAL::collinear(
-                                query_point, adj_points[i], closest_point));
-                    break;
-                }
-            }
+      }
+    }
+    if (d == std::numeric_limits<size_t>::max()) {
+      Eigen::MatrixXd tmp_vertices(V.rows(), V.cols());
+      for (size_t i=0; i<V.rows(); i++) {
+        for (size_t j=0; j<V.cols(); j++) {
+          tmp_vertices(i,j) = CGAL::to_double(V(i,j));
         }
-        assert(d != std::numeric_limits<size_t>::max());
-
-        return process_edge_case(query_idx, s, d, preferred_facet, orientation);
-    };
-
-    const size_t num_queries = P.rows();
-    R.resize(num_queries, 1);
-    S.resize(num_queries, 1);
-    for (size_t i=0; i<num_queries; i++) {
-        const Point_3 query(P(i,0), P(i,1), P(i,2));
-        auto projection = tree.closest_point_and_primitive(query);
-        const Point_3 closest_point = projection.first;
-        size_t fid = projection.second - triangles.begin();
-        bool fid_ori = false;
-
-        // Gether all facets sharing the closest point.
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(closest_point, query),
-                std::back_inserter(intersected_faces));
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-
-        size_t element_index;
-        auto element_type = determine_element_type(closest_point, fid,
-                element_index);
-        switch(element_type) {
-            case VERTEX:
-                {
-                    const auto& f = F.row(I(fid, 0));
-                    const size_t s = f[element_index];
-                    fid = process_vertex_case(i, s, I(fid, 0), fid_ori);
-                }
-                break;
-            case EDGE:
-                {
-                    const auto& f = F.row(I(fid, 0));
-                    const size_t s = f[(element_index+1)%3];
-                    const size_t d = f[(element_index+2)%3];
-                    fid = process_edge_case(i, s, d, I(fid, 0), fid_ori);
-                }
-                break;
-            case FACE:
-                {
-                    fid = process_face_case(i, fid, fid_ori);
-                }
-                break;
-            default:
-                throw std::runtime_error("Unknown element type.");
+      }
+      Eigen::MatrixXi tmp_faces(adj_faces.size(), 3);
+      for (size_t i=0; i<adj_faces.size(); i++) {
+        tmp_faces.row(i) = F.row(adj_faces[i]);
+      }
+      igl::writePLY("debug.ply", tmp_vertices, tmp_faces, false);
+      throw std::runtime_error("Invalid vertex neighborhood");
+    }
+    const auto itr = v2f.equal_range(d);
+    assert(itr.first != itr.second);
+
+    return process_edge_case(query_idx, s, d, itr.first->second, orientation);
+  };
+
+  const size_t num_queries = P.rows();
+  R.resize(num_queries, 1);
+  S.resize(num_queries, 1);
+  for (size_t i=0; i<num_queries; i++) {
+    const Point_3 query(P(i,0), P(i,1), P(i,2));
+    auto projection = tree.closest_point_and_primitive(query);
+    const Point_3 closest_point = projection.first;
+    size_t fid = projection.second - triangles.begin();
+    bool fid_ori = false;
+
+    // Gether all facets sharing the closest point.
+    std::vector<Tree::Primitive_id> intersected_faces;
+    tree.all_intersected_primitives(Segment_3(closest_point, query),
+        std::back_inserter(intersected_faces));
+    const size_t num_intersected_faces = intersected_faces.size();
+    std::vector<size_t> intersected_face_indices(num_intersected_faces);
+    std::transform(intersected_faces.begin(),
+        intersected_faces.end(),
+        intersected_face_indices.begin(),
+        [&](const Tree::Primitive_id& itr) -> size_t
+        { return I(itr-triangles.begin(), 0); });
+
+    size_t element_index;
+    auto element_type = determine_element_type(closest_point, fid,
+        element_index);
+    switch(element_type) {
+      case VERTEX:
+        {
+          const auto& f = F.row(I(fid, 0));
+          const size_t s = f[element_index];
+          fid = process_vertex_case(i, s, I(fid, 0), fid_ori);
+        }
+        break;
+      case EDGE:
+        {
+          const auto& f = F.row(I(fid, 0));
+          const size_t s = f[(element_index+1)%3];
+          const size_t d = f[(element_index+2)%3];
+          fid = process_edge_case(i, s, d, I(fid, 0), fid_ori);
         }
+        break;
+      case FACE:
+        {
+          fid = process_face_case(i, closest_point, fid, fid_ori);
+        }
+        break;
+      default:
+        throw std::runtime_error("Unknown element type.");
+    }
 
 
-        R(i,0) = fid;
-        S(i,0) = fid_ori;
-    }
+    R(i,0) = fid;
+    S(i,0) = fid_ori;
+  }
 }
 
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedP,
-    typename DerivedR,
-    typename DerivedS >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedR,
+  typename DerivedS >
 IGL_INLINE void igl::copyleft::cgal::closest_facet(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        Eigen::PlainObjectBase<DerivedR>& R,
-        Eigen::PlainObjectBase<DerivedS>& S) {
-    const size_t num_faces = F.rows();
-    Eigen::VectorXi I(num_faces);
-    I.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::copyleft::cgal::closest_facet(V, F, I, P, R, S);
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S) {
+  const size_t num_faces = F.rows();
+  Eigen::VectorXi I(num_faces);
+  I.setLinSpaced(num_faces, 0, num_faces-1);
+  igl::copyleft::cgal::closest_facet(V, F, I, P, uE2E, EMAP, R, S);
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::copyleft::cgal::closest_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::closest_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/closest_facet.h

@@ -11,6 +11,7 @@
 
 #include "../../igl_inline.h"
 #include <Eigen/Core>
+#include <vector>
 
 namespace igl
 {
@@ -35,6 +36,8 @@ namespace igl
           typename DerivedF,
           typename DerivedI,
           typename DerivedP,
+          typename uE2EType,
+          typename DerivedEMAP,
           typename DerivedR,
           typename DerivedS >
       IGL_INLINE void closest_facet(
@@ -42,6 +45,8 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedI>& I,
         const Eigen::PlainObjectBase<DerivedP>& P,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
 
@@ -49,12 +54,16 @@ namespace igl
           typename DerivedV,
           typename DerivedF,
           typename DerivedP,
+          typename uE2EType,
+          typename DerivedEMAP,
           typename DerivedR,
           typename DerivedS >
       IGL_INLINE void closest_facet(
               const Eigen::PlainObjectBase<DerivedV>& V,
               const Eigen::PlainObjectBase<DerivedF>& F,
               const Eigen::PlainObjectBase<DerivedP>& P,
+              const std::vector<std::vector<uE2EType> >& uE2E,
+              const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
     }

include/igl/copyleft/cgal/extract_cells.cpp

@@ -229,7 +229,7 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
             const auto& I = Is[i];
             Eigen::VectorXi closest_facets, closest_facet_orientations;
             igl::copyleft::cgal::closest_facet(V, F, I, queries,
-                    closest_facets, closest_facet_orientations);
+                    uE2E, EMAP, closest_facets, closest_facet_orientations);
 
             for (size_t j=0; j<num_components; j++) {
                 if (i == j) continue;

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -36,30 +36,6 @@ namespace igl {
             typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
             typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
 
-            enum ElementType { VERTEX, EDGE, FACE };
-            template<typename DerivedV, typename DerivedF, typename DerivedI>
-            ElementType determine_element_type(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedI>& I,
-                    const size_t fid, const Point_3& p,
-                    size_t& element_index) {
-                const Eigen::Vector3i f = F.row(I(fid, 0));
-                const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
-                const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
-                const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
-
-                if (p == p0) { element_index = 0; return VERTEX; }
-                if (p == p1) { element_index = 1; return VERTEX; }
-                if (p == p2) { element_index = 2; return VERTEX; }
-                if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
-                if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
-                if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
-
-                element_index = 0;
-                return FACE;
-            }
-
             template<typename DerivedF, typename DerivedI>
             void extract_adj_faces(
                     const Eigen::PlainObjectBase<DerivedF>& F,
@@ -296,6 +272,25 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
     Tree tree(triangles.begin(), triangles.end());
     tree.accelerate_distance_queries();
 
+    enum ElementType { VERTEX, EDGE, FACE };
+    auto determine_element_type = [&](
+            size_t fid, const Point_3& p, size_t& element_index) -> ElementType{
+        const Eigen::Vector3i f = F.row(I(fid, 0));
+        const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+        const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+        const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+
+        if (p == p0) { element_index = 0; return VERTEX; }
+        if (p == p1) { element_index = 1; return VERTEX; }
+        if (p == p2) { element_index = 2; return VERTEX; }
+        if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+        if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+        if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+        element_index = 0;
+        return FACE;
+    };
+
     const size_t num_queries = P.rows();
     inside.resize(num_queries, 1);
     for (size_t i=0; i<num_queries; i++) {
@@ -305,8 +300,7 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
         size_t fid = projection.second - triangles.begin();
 
         size_t element_index;
-        switch (determine_element_type(
-                    V, F, I, fid, closest_point, element_index)) {
+        switch (determine_element_type(fid, closest_point, element_index)) {
             case VERTEX:
                 {
                     const size_t s = F(I(fid, 0), element_index);

include/igl/copyleft/cgal/points_inside_component.h

@@ -17,41 +17,49 @@ namespace igl {
   {
     namespace cgal {
 
-        // Determine if queries points P are inside of connected facet component
-        // (V, F, I), where I indicates a subset of facets that forms the
-        // component.
-        //
-        // Precondition:
-        // The input mesh must be a closed, self-intersection free,
-        // non-degenerated surface.  Queries points must be either inside or
-        // outside of the mesh (i.e. not on the surface of the mesh).
-        //
-        // Inputs:
-        //   V  #V by 3 array of vertex positions.
-        //   F  #F by 3 array of triangles.
-        //   I  #I list of triangle indices to consider.
-        //   P  #P by 3 array of query points.
-        //
-        // Outputs:
-        //   inside  #P list of booleans that is true iff the corresponding
-        //           query point is inside of the mesh.
-        template<typename DerivedV, typename DerivedF, typename DerivedI,
-            typename DerivedP, typename DerivedB>
-            IGL_INLINE void points_inside_component(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedI>& I,
-                    const Eigen::PlainObjectBase<DerivedP>& P,
-                    Eigen::PlainObjectBase<DerivedB>& inside);
+      // Determine if queries points P are inside of connected facet component
+      // (V, F, I), where I indicates a subset of facets that forms the
+      // component.
+      //
+      // Precondition:
+      // The input mesh must be a closed, self-intersection free,
+      // non-degenerated surface.  Queries points must be either inside or
+      // outside of the mesh (i.e. not on the surface of the mesh).
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertex positions.
+      //   F  #F by 3 array of triangles.
+      //   I  #I list of triangle indices to consider.
+      //   P  #P by 3 array of query points.
+      //
+      // Outputs:
+      //   inside  #P list of booleans that is true iff the corresponding
+      //           query point is inside of the mesh.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename DerivedP,
+        typename DerivedB>
+      IGL_INLINE void points_inside_component(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside);
 
-        // Determine if query points P is inside of the mesh (V, F).
-        // See above for precondition and I/O specs.
-        template<typename DerivedV, typename DerivedF, typename DerivedP, typename DerivedB>
-            IGL_INLINE void points_inside_component(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedP>& P,
-                    Eigen::PlainObjectBase<DerivedB>& inside);
+      // Determine if query points P is inside of the mesh (V, F).
+      // See above for precondition and I/O specs.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename DerivedB>
+      IGL_INLINE void points_inside_component(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          Eigen::PlainObjectBase<DerivedB>& inside);
     }
   }
 }

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -26,235 +26,234 @@
 #include <queue>
 
 namespace propagate_winding_numbers_helper {
-    template<
-        typename DerivedF,
-        typename DeriveduE,
-        typename uE2EType >
-    bool is_orientable(
-            const Eigen::PlainObjectBase<DerivedF>& F,
-            const Eigen::PlainObjectBase<DeriveduE>& uE,
-            const std::vector<std::vector<uE2EType> >& uE2E) {
-        const size_t num_faces = F.rows();
-        const size_t num_edges = uE.rows();
-        auto edge_index_to_face_index = [&](size_t ei) {
-            return ei % num_faces;
-        };
-        auto is_consistent = [&](size_t fid, size_t s, size_t d) {
-            if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return true;
-            if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return true;
-            if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return true;
+  template<
+    typename DerivedF,
+    typename DeriveduE,
+    typename uE2EType >
+  bool is_orientable(
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      const Eigen::PlainObjectBase<DeriveduE>& uE,
+      const std::vector<std::vector<uE2EType> >& uE2E) {
+    const size_t num_faces = F.rows();
+    const size_t num_edges = uE.rows();
+    auto edge_index_to_face_index = [&](size_t ei) {
+      return ei % num_faces;
+    };
+    auto is_consistent = [&](size_t fid, size_t s, size_t d) {
+      if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return true;
+      if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return true;
+      if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return true;
 
-            if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return false;
-            if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return false;
-            if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return false;
-            throw "Invalid face!!";
-        };
-        for (size_t i=0; i<num_edges; i++) {
-            const size_t s = uE(i,0);
-            const size_t d = uE(i,1);
-            int count=0;
-            for (const auto& ei : uE2E[i]) {
-                const size_t fid = edge_index_to_face_index(ei);
-                if (is_consistent(fid, s, d)) count++;
-                else count--;
-            }
-            if (count != 0) {
-                return false;
-            }
-        }
-        return true;
+      if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return false;
+      if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return false;
+      if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return false;
+      throw "Invalid face!!";
+    };
+    for (size_t i=0; i<num_edges; i++) {
+      const size_t s = uE(i,0);
+      const size_t d = uE(i,1);
+      int count=0;
+      for (const auto& ei : uE2E[i]) {
+        const size_t fid = edge_index_to_face_index(ei);
+        if (is_consistent(fid, s, d)) count++;
+        else count--;
+      }
+      if (count != 0) {
+        return false;
+      }
     }
+    return true;
+  }
 }
 
 template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedL,
-typename DerivedW>
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedL,
+  typename DerivedW>
 IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedL>& labels,
-        Eigen::PlainObjectBase<DerivedW>& W) {
-    const size_t num_faces = F.rows();
-    //typedef typename DerivedF::Scalar Index;
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedL>& labels,
+    Eigen::PlainObjectBase<DerivedW>& W) {
+  const size_t num_faces = F.rows();
+  //typedef typename DerivedF::Scalar Index;
+
+  Eigen::MatrixXi E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  if (!propagate_winding_numbers_helper::is_orientable(F, uE, uE2E)) {
+      std::cerr << "Input mesh is not orientable!" << std::endl;
+  }
 
-    Eigen::MatrixXi E, uE;
-    Eigen::VectorXi EMAP;
-    std::vector<std::vector<size_t> > uE2E;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-    assert(propagate_winding_numbers_helper::is_orientable(F, uE, uE2E));
+  Eigen::VectorXi P;
+  const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
 
-    Eigen::VectorXi P;
-    const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+  DerivedW per_patch_cells;
+  const size_t num_cells =
+    igl::copyleft::cgal::extract_cells(
+        V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
 
-    DerivedW per_patch_cells;
-    const size_t num_cells =
-        igl::copyleft::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+  typedef std::tuple<size_t, bool, size_t> CellConnection;
+  std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+  for (size_t i=0; i<num_patches; i++) {
+    const int positive_cell = per_patch_cells(i,0);
+    const int negative_cell = per_patch_cells(i,1);
+    cell_adjacency[positive_cell].emplace(negative_cell, false, i);
+    cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+  }
 
-    typedef std::tuple<size_t, bool, size_t> CellConnection;
-    std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+  auto save_cell = [&](const std::string& filename, size_t cell_id) {
+    std::vector<size_t> faces;
     for (size_t i=0; i<num_patches; i++) {
-        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);
+      if ((per_patch_cells.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);
+          }
+        }
+      }
+    }
+    Eigen::MatrixXi cell_faces(faces.size(), 3);
+    for (size_t i=0; i<faces.size(); i++) {
+      cell_faces.row(i) = F.row(faces[i]);
     }
+    Eigen::MatrixXd vertices(V.rows(), 3);
+    for (size_t i=0; i<(size_t)V.rows(); i++) {
+      assign_scalar(V(i,0), vertices(i,0));
+      assign_scalar(V(i,1), vertices(i,1));
+      assign_scalar(V(i,2), vertices(i,2));
+    }
+    writePLY(filename, vertices, cell_faces);
+  };
 
-    auto save_cell = [&](const std::string& filename, size_t cell_id) {
-        std::vector<size_t> faces;
-        for (size_t i=0; i<num_patches; i++) {
-            if ((per_patch_cells.row(i).array() == cell_id).any()) {
-                for (size_t j=0; j<num_faces; j++) {
-                    if ((size_t)P[j] == i) {
-                        faces.push_back(j);
-                    }
-                }
-            }
-        }
-        Eigen::MatrixXi cell_faces(faces.size(), 3);
-        for (size_t i=0; i<faces.size(); i++) {
-            cell_faces.row(i) = F.row(faces[i]);
-        }
-        Eigen::MatrixXd vertices(V.rows(), 3);
-        for (size_t i=0; i<(size_t)V.rows(); i++) {
-            assign_scalar(V(i,0), vertices(i,0));
-            assign_scalar(V(i,1), vertices(i,1));
-            assign_scalar(V(i,2), vertices(i,2));
-        }
-        writePLY(filename, vertices, cell_faces);
+#ifndef NDEBUG
+  {
+    // Check for odd cycle.
+    Eigen::VectorXi cell_labels(num_cells);
+    cell_labels.setZero();
+    Eigen::VectorXi parents(num_cells);
+    parents.setConstant(-1);
+    auto trace_parents = [&](size_t idx) {
+      std::list<size_t> path;
+      path.push_back(idx);
+      while ((size_t)parents[path.back()] != path.back()) {
+        path.push_back(parents[path.back()]);
+      }
+      return path;
     };
-
-    {
-        // Check for odd cycle.
-        Eigen::VectorXi cell_labels(num_cells);
-        cell_labels.setZero();
-        Eigen::VectorXi parents(num_cells);
-        parents.setConstant(-1);
-        auto trace_parents = [&](size_t idx) {
-            std::list<size_t> path;
-            path.push_back(idx);
-            while ((size_t)parents[path.back()] != path.back()) {
-                path.push_back(parents[path.back()]);
-            }
-            return path;
-        };
-        for (size_t i=0; i<num_cells; i++) {
-            if (cell_labels[i] == 0) {
-                cell_labels[i] = 1;
-                std::queue<size_t> Q;
-                Q.push(i);
-                parents[i] = i;
-                while (!Q.empty()) {
-                    size_t curr_idx = Q.front();
-                    Q.pop();
-                    int curr_label = cell_labels[curr_idx];
-                    for (const auto& neighbor : cell_adjacency[curr_idx]) {
-                        if (cell_labels[std::get<0>(neighbor)] == 0) {
-                            cell_labels[std::get<0>(neighbor)] = curr_label * -1;
-                            Q.push(std::get<0>(neighbor));
-                            parents[std::get<0>(neighbor)] = curr_idx;
-                        } else {
-                            if (cell_labels[std::get<0>(neighbor)] !=
-                                    curr_label * -1) {
-                                std::cerr << "Odd cell cycle detected!" << std::endl;
-                                auto path = trace_parents(curr_idx);
-                                path.reverse();
-                                auto path2 = trace_parents(std::get<0>(neighbor));
-                                path.insert(path.end(),
-                                        path2.begin(), path2.end());
-                                for (auto cell_id : path) {
-                                    std::cout << cell_id << " ";
-                                    std::stringstream filename;
-                                    filename << "cell_" << cell_id << ".ply";
-                                    save_cell(filename.str(), cell_id);
-                                }
-                                std::cout << std::endl;
-                            }
-                            assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
-                        }
-                    }
+    for (size_t i=0; i<num_cells; i++) {
+      if (cell_labels[i] == 0) {
+        cell_labels[i] = 1;
+        std::queue<size_t> Q;
+        Q.push(i);
+        parents[i] = i;
+        while (!Q.empty()) {
+          size_t curr_idx = Q.front();
+          Q.pop();
+          int curr_label = cell_labels[curr_idx];
+          for (const auto& neighbor : cell_adjacency[curr_idx]) {
+            if (cell_labels[std::get<0>(neighbor)] == 0) {
+              cell_labels[std::get<0>(neighbor)] = curr_label * -1;
+              Q.push(std::get<0>(neighbor));
+              parents[std::get<0>(neighbor)] = curr_idx;
+            } else {
+              if (cell_labels[std::get<0>(neighbor)] !=
+                  curr_label * -1) {
+                std::cerr << "Odd cell cycle detected!" << std::endl;
+                auto path = trace_parents(curr_idx);
+                path.reverse();
+                auto path2 = trace_parents(std::get<0>(neighbor));
+                path.insert(path.end(),
+                    path2.begin(), path2.end());
+                for (auto cell_id : path) {
+                  std::cout << cell_id << " ";
+                  std::stringstream filename;
+                  filename << "cell_" << cell_id << ".ply";
+                  save_cell(filename.str(), cell_id);
                 }
+                std::cout << std::endl;
+              }
+              assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
             }
+          }
         }
+      }
     }
+  }
+#endif
 
-    size_t outer_facet;
-    bool flipped;
-    Eigen::VectorXi I;
-    I.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
+  size_t outer_facet;
+  bool flipped;
+  Eigen::VectorXi I;
+  I.setLinSpaced(num_faces, 0, num_faces-1);
+  igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
 
-    const size_t outer_patch = P[outer_facet];
-    const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
+  const size_t outer_patch = P[outer_facet];
+  const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
 
-    Eigen::VectorXi patch_labels(num_patches);
-    const int INVALID = std::numeric_limits<int>::max();
-    patch_labels.setConstant(INVALID);
-    for (size_t i=0; i<num_faces; i++) {
-        if (patch_labels[P[i]] == INVALID) {
-            patch_labels[P[i]] = labels[i];
-        } else {
-            assert(patch_labels[P[i]] == labels[i]);
-        }
+  Eigen::VectorXi patch_labels(num_patches);
+  const int INVALID = std::numeric_limits<int>::max();
+  patch_labels.setConstant(INVALID);
+  for (size_t i=0; i<num_faces; i++) {
+    if (patch_labels[P[i]] == INVALID) {
+      patch_labels[P[i]] = labels[i];
+    } else {
+      assert(patch_labels[P[i]] == labels[i]);
     }
-    assert((patch_labels.array() != INVALID).all());
-    const size_t num_labels = patch_labels.maxCoeff()+1;
+  }
+  assert((patch_labels.array() != INVALID).all());
+  const size_t num_labels = patch_labels.maxCoeff()+1;
 
-    Eigen::MatrixXi per_cell_W(num_cells, num_labels);
-    per_cell_W.setConstant(INVALID);
-    per_cell_W.row(infinity_cell).setZero();
-    std::queue<size_t> Q;
-    Q.push(infinity_cell);
-    while (!Q.empty()) {
-        size_t curr_cell = Q.front();
-        Q.pop();
-        for (const auto& neighbor : cell_adjacency[curr_cell]) {
-            size_t neighbor_cell, patch_idx;
-            bool direction;
-            std::tie(neighbor_cell, direction, patch_idx) = neighbor;
-            if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
-                per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
-                for (size_t i=0; i<num_labels; i++) {
-                    int inc = (patch_labels[patch_idx] == (int)i) ?
-                        (direction ? -1:1) :0;
-                    per_cell_W(neighbor_cell, i) =
-                        per_cell_W(curr_cell, i) + inc;
-                }
-                Q.push(neighbor_cell);
-            } else {
-                for (size_t i=0; i<num_labels; i++) {
-                    if ((int)i == patch_labels[patch_idx]) {
+  Eigen::MatrixXi per_cell_W(num_cells, num_labels);
+  per_cell_W.setConstant(INVALID);
+  per_cell_W.row(infinity_cell).setZero();
+  std::queue<size_t> Q;
+  Q.push(infinity_cell);
+  while (!Q.empty()) {
+    size_t curr_cell = Q.front();
+    Q.pop();
+    for (const auto& neighbor : cell_adjacency[curr_cell]) {
+      size_t neighbor_cell, patch_idx;
+      bool direction;
+      std::tie(neighbor_cell, direction, patch_idx) = neighbor;
+      if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
+        per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
+        for (size_t i=0; i<num_labels; i++) {
+          int inc = (patch_labels[patch_idx] == (int)i) ?
+            (direction ? -1:1) :0;
+          per_cell_W(neighbor_cell, i) =
+            per_cell_W(curr_cell, i) + inc;
+        }
+        Q.push(neighbor_cell);
+      } else {
 #ifndef NDEBUG
-                        int inc = direction ? -1:1;
-                        assert(per_cell_W(neighbor_cell, i) ==
-                                per_cell_W(curr_cell, i) + inc);
-#endif
-                    } else {
-                        assert(per_cell_W(neighbor_cell, i) ==
-                                per_cell_W(curr_cell, i));
-                    }
-                }
-            }
+        for (size_t i=0; i<num_labels; i++) {
+          if ((int)i == patch_labels[patch_idx]) {
+            int inc = direction ? -1:1;
+            assert(per_cell_W(neighbor_cell, i) ==
+                per_cell_W(curr_cell, i) + inc);
+          } else {
+            assert(per_cell_W(neighbor_cell, i) ==
+                per_cell_W(curr_cell, i));
+          }
         }
+#endif
+      }
     }
+  }
 
-    W.resize(num_faces, num_labels*2);
-    for (size_t i=0; i<num_faces; i++) {
-        const size_t patch = P[i];
-        const size_t positive_cell = per_patch_cells(patch, 0);
-        const size_t negative_cell = per_patch_cells(patch, 1);
-        for (size_t j=0; j<num_labels; j++) {
-            W(i,j*2  ) = per_cell_W(positive_cell, j);
-            W(i,j*2+1) = per_cell_W(negative_cell, j);
-        }
+  W.resize(num_faces, num_labels*2);
+  for (size_t i=0; i<num_faces; i++) {
+    const size_t patch = P[i];
+    const size_t positive_cell = per_patch_cells(patch, 0);
+    const size_t negative_cell = per_patch_cells(patch, 1);
+    for (size_t j=0; j<num_labels; j++) {
+      W(i,j*2  ) = per_cell_W(positive_cell, j);
+      W(i,j*2+1) = per_cell_W(negative_cell, j);
     }
-
-    //for (size_t i=0; i<num_cells; i++) {
-    //    std::stringstream filename;
-    //    filename << "cell_" << i << "_w_" << per_cell_W(i, 1) << ".ply";
-    //    save_cell(filename.str(), i);
-    //}
+  }
 }
 
 

include/igl/copyleft/cgal/propagate_winding_numbers.h

@@ -24,32 +24,32 @@ namespace igl {
   {
     namespace cgal {
 
-        // Compute winding number on each side of the face.  The input mesh
-        // could contain multiple connected components.  The input mesh must
-        // represent the boundary of a valid 3D volume, which means it is
-        // closed, consistently oriented and induces integer winding numbers.
-        //
-        // Inputs:
-        //   V  #V by 3 list of vertex positions.
-        //   F  #F by 3 list of triangle indices into V.
-        //   labels  #F list of facet labels ranging from 0 to k-1.
-        //
-        // Output:
-        //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
-        //      number on the positive side of facet ``i`` with respect to the
-        //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
-        //      number on the negative side of facet ``i`` with respect to the
-        //      facets labeled ``j``.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedL,
-            typename DerivedW>
-        IGL_INLINE void propagate_winding_numbers(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedL>& labels,
-                Eigen::PlainObjectBase<DerivedW>& W);
+      // Compute winding number on each side of the face.  The input mesh
+      // could contain multiple connected components.  The input mesh must
+      // represent the boundary of a valid 3D volume, which means it is
+      // closed, consistently oriented and induces integer winding numbers.
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions.
+      //   F  #F by 3 list of triangle indices into V.
+      //   labels  #F list of facet labels ranging from 0 to k-1.
+      //
+      // Output:
+      //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
+      //      number on the positive side of facet ``i`` with respect to the
+      //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
+      //      number on the negative side of facet ``i`` with respect to the
+      //      facets labeled ``j``.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedL,
+        typename DerivedW>
+      IGL_INLINE void propagate_winding_numbers(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const Eigen::PlainObjectBase<DerivedL>& labels,
+          Eigen::PlainObjectBase<DerivedW>& W);
     }
   }
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -272,7 +272,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
 
     // Process un-touched faces.
     for (size_t i=0; i<num_faces; i++) {
-        if (!is_offending[i]) {
+        if (!is_offending[i] && !T[i].is_degenerate()) {
             resolved_faces.push_back(
                     { F(i,0), F(i,1), F(i,2) } );
             source_faces.push_back(i);

include/igl/vertex_triangle_adjacency.cpp

@@ -50,4 +50,5 @@ template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 1, -1,
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, unsigned long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > >&, std::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > > >&);
 #endif