Merge pull request #126 from qnzhou/master

Robust inside/outside test

Former-commit-id: 7c9366dbae7b4a4ecf46ad6bb5240415c8c15ec9

include/igl/cgal/is_inside.cpp

@@ -0,0 +1,367 @@
+#include "is_inside.h"
+
+#include <cassert>
+#include <list>
+#include <limits>
+#include <vector>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include "order_facets_around_edge.h"
+#include "assign_scalar.h"
+
+namespace igl {
+    namespace cgal {
+        namespace is_inside_helper {
+            typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+            typedef Kernel::Ray_3 Ray_3;
+            typedef Kernel::Point_3 Point_3;
+            typedef Kernel::Vector_3 Vector_3;
+            typedef Kernel::Triangle_3 Triangle;
+            typedef Kernel::Plane_3 Plane_3;
+            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;
+
+            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 DerivedV, typename DerivedF, typename DerivedI>
+            void extract_adj_faces(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const size_t s, const size_t d,
+                    std::vector<int>& adj_faces) {
+                const size_t num_faces = I.rows();
+                for (size_t i=0; i<num_faces; i++) {
+                    Eigen::Vector3i f = F.row(I(i, 0));
+                    if ((f[0] == s && f[1] == d) ||
+                        (f[1] == s && f[2] == d) ||
+                        (f[2] == s && f[0] == d)) {
+                        adj_faces.push_back((I(i, 0)+1) * -1);
+                        continue;
+                    }
+                    if ((f[0] == d && f[1] == s) ||
+                        (f[1] == d && f[2] == s) ||
+                        (f[2] == d && f[0] == s)) {
+                        adj_faces.push_back(I(i, 0)+1);
+                        continue;
+                    }
+                }
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            void extract_adj_vertices(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const size_t v, std::vector<int>& adj_vertices) {
+                std::set<size_t> unique_adj_vertices;
+                const size_t num_faces = I.rows();
+                for (size_t i=0; i<num_faces; i++) {
+                    Eigen::Vector3i f = F.row(I(i, 0));
+                    assert((f.array() < V.rows()).all());
+                    if (f[0] == v) {
+                        unique_adj_vertices.insert(f[1]);
+                        unique_adj_vertices.insert(f[2]);
+                    } else if (f[1] == v) {
+                        unique_adj_vertices.insert(f[0]);
+                        unique_adj_vertices.insert(f[2]);
+                    } else if (f[2] == v) {
+                        unique_adj_vertices.insert(f[0]);
+                        unique_adj_vertices.insert(f[1]);
+                    }
+                }
+                adj_vertices.resize(unique_adj_vertices.size());
+                std::copy(unique_adj_vertices.begin(),
+                        unique_adj_vertices.end(),
+                        adj_vertices.begin());
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_edge_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t s, size_t d) {
+                // Algorithm:
+                //
+                // Order the adj faces around the edge (s,d) clockwise using
+                // query point as pivot.  (i.e. The first face of the ordering
+                // is directly after the pivot point, and the last face is
+                // directly before the pivot.)
+                //
+                // The point is outside if the first and last faces of the
+                // ordering forms a convex angle.  This check can be done
+                // without any construction by looking at the orientation of the
+                // faces.  The angle is convex iff the first face contains (s,d)
+                // as an edge and the last face contains (d,s) as an edge.
+                //
+                // The point is inside if the first and last faces of the
+                // ordering forms a concave angle.  That is the first face
+                // contains (d,s) as an edge and the last face contains (s,d) as
+                // an edge.
+                //
+                // In the special case of duplicated faces. I.e. multiple faces
+                // sharing the same 3 corners, but not necessarily the same
+                // orientation.  The ordering will always rank faces containing
+                // edge (s,d) before faces containing edge (d,s).
+                //
+                // Therefore, if there are any duplicates of the first faces,
+                // the ordering will always choose the one with edge (s,d) if
+                // possible.  The same for the last face.
+                //
+                // In the very degenerated case where the first and last face
+                // are duplicates, but with different orientations, it is
+                // equally valid to think the angle formed by them is either 0
+                // or 360 degrees.  By default, 0 degree is used, and thus the
+                // query point is outside.
+
+                std::vector<int> adj_faces;
+                extract_adj_faces(V, F, I, s, d, adj_faces);
+                const size_t num_adj_faces = adj_faces.size();
+                assert(num_adj_faces > 0);
+
+                DerivedV pivot_point(1, 3);
+                igl::cgal::assign_scalar(query.x(), pivot_point(0, 0));
+                igl::cgal::assign_scalar(query.y(), pivot_point(0, 1));
+                igl::cgal::assign_scalar(query.z(), pivot_point(0, 2));
+                Eigen::VectorXi order;
+                order_facets_around_edge(V, F, s, d,
+                        adj_faces, pivot_point, order);
+                assert(order.size() == num_adj_faces);
+                if (adj_faces[order[0]] > 0 &&
+                    adj_faces[order[num_adj_faces-1] < 0]) {
+                    return true;
+                } else if (adj_faces[order[0]] < 0 &&
+                    adj_faces[order[num_adj_faces-1] > 0]) {
+                    return false;
+                } else {
+                    throw "The input mesh does not represent a valid volume";
+                }
+                assert(false);
+                return false;
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_vertex_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t s) {
+                std::vector<int> adj_vertices;
+                extract_adj_vertices(V, F, I, s, adj_vertices);
+                const size_t num_adj_vertices = adj_vertices.size();
+
+                std::vector<Point_3> adj_points;
+                for (size_t i=0; i<num_adj_vertices; i++) {
+                    const size_t vi = adj_vertices[i];
+                    adj_points.emplace_back(V(vi,0), V(vi,1), V(vi,2));
+                }
+
+                // A plane is on the exterior if all adj_points lies on or to
+                // one side of the plane.
+                auto is_on_exterior = [&](const Plane_3& separator) {
+                    size_t positive=0;
+                    size_t negative=0;
+                    size_t coplanar=0;
+                    for (const auto& point : adj_points) {
+                        switch(separator.oriented_side(point)) {
+                            case CGAL::ON_POSITIVE_SIDE:
+                                positive++;
+                                break;
+                            case CGAL::ON_NEGATIVE_SIDE:
+                                negative++;
+                                break;
+                            case CGAL::ON_ORIENTED_BOUNDARY:
+                                coplanar++;
+                                break;
+                            default:
+                                assert(false);
+                        }
+                    }
+                    auto query_orientation = separator.oriented_side(query);
+                    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();
+                Point_3 p(V(s,0), V(s,1), V(s,2));
+                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++) {
+                        const size_t vj = adj_vertices[j];
+                        Plane_3 separator(p, adj_points[i], adj_points[j]);
+                        assert(!separator.is_degenerate());
+                        if (is_on_exterior(separator)) {
+                            d = vi;
+                            assert(!CGAL::collinear(p, adj_points[i], query));
+                            break;
+                        }
+                    }
+                    if (d < V.rows()) break;
+                }
+                if (d > V.rows()) {
+                    // All adj faces are coplanar, use the first edge.
+                    d = adj_vertices[0];
+                }
+                return determine_point_edge_orientation(V, F, I, query, s, d);
+            }
+
+            template<typename DerivedV, typename DerivedF, typename DerivedI>
+            bool determine_point_face_orientation(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Point_3& query, size_t fid) {
+                // Algorithm: A point is on the inside of a face if the
+                // tetrahedron formed by them is negatively oriented.
+
+                Eigen::Vector3i f = F.row(I(fid, 0));
+                const Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+                const Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+                const Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+                auto result = CGAL::orientation(v0, v1, v2, query);
+                assert(result != CGAL::COPLANAR);
+                return result == CGAL::NEGATIVE;
+            }
+        }
+    }
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedI>
+IGL_INLINE bool igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V1,
+        const Eigen::PlainObjectBase<DerivedF>& F1,
+        const Eigen::PlainObjectBase<DerivedI>& I1,
+        const Eigen::PlainObjectBase<DerivedV>& V2,
+        const Eigen::PlainObjectBase<DerivedF>& F2,
+        const Eigen::PlainObjectBase<DerivedI>& I2) {
+    assert(F1.rows() > 0);
+    assert(I1.rows() > 0);
+    assert(F2.rows() > 0);
+    assert(I2.rows() > 0);
+
+    const Eigen::Vector3i& f = F1.row(I1(0, 0));
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> query(
+            (V1(f[0],0) + V1(f[1],0) + V1(f[2],0))/3.0,
+            (V1(f[0],1) + V1(f[1],1) + V1(f[2],1))/3.0,
+            (V1(f[0],2) + V1(f[1],2) + V1(f[2],2))/3.0);
+    Eigen::VectorXi inside;
+    igl::cgal::is_inside_batch(V2, F2, I2, query, inside);
+    assert(inside.size() == 1);
+    return inside[0];
+}
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::cgal::is_inside(
+        const Eigen::PlainObjectBase<DerivedV>& V1,
+        const Eigen::PlainObjectBase<DerivedF>& F1,
+        const Eigen::PlainObjectBase<DerivedV>& V2,
+        const Eigen::PlainObjectBase<DerivedF>& F2) {
+    Eigen::VectorXi I1(F1.rows()), I2(F2.rows());
+    I1.setLinSpaced(F1.rows(), 0, F1.rows()-1);
+    I2.setLinSpaced(F2.rows(), 0, F2.rows()-1);
+    return igl::cgal::is_inside(V1, F1, I1, V2, F2, I2);
+}
+
+template<typename DerivedV, typename DerivedF, typename DerivedI,
+    typename DerivedP, typename DerivedB>
+IGL_INLINE void igl::cgal::is_inside_batch(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside) {
+    using namespace igl::cgal::is_inside_helper;
+    assert(F.rows() > 0);
+    assert(I.rows() > 0);
+
+    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)));
+        assert(!triangles.back().is_degenerate());
+    }
+    Tree tree(triangles.begin(), triangles.end());
+    tree.accelerate_distance_queries();
+
+    const size_t num_queries = P.rows();
+    inside.resize(num_queries);
+    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);
+        auto closest_point = projection.first;
+        size_t fid = projection.second - triangles.begin();
+
+        size_t element_index;
+        switch (determine_element_type(
+                    V, F, I, fid, closest_point, element_index)) {
+            case VERTEX:
+                {
+                    const size_t s = F(I(fid, 0), element_index);
+                    inside[i] = determine_point_vertex_orientation(
+                            V, F, I, query, s);
+                }
+                break;
+            case EDGE:
+                {
+                    const size_t s = F(I(fid, 0), (element_index+1)%3);
+                    const size_t d = F(I(fid, 0), (element_index+2)%3);
+                    inside[i] = determine_point_edge_orientation(
+                            V, F, I, query, s, d);
+                }
+                break;
+            case FACE:
+                inside[i] = determine_point_face_orientation(V, F, I, query, fid);
+                break;
+            default:
+                assert(false);
+        }
+    }
+}
+
+template<typename DerivedV, typename DerivedF, typename DerivedP,
+    typename DerivedB>
+IGL_INLINE void igl::cgal::is_inside_batch(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside) {
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+    igl::cgal::is_inside_batch(V, F, I, P, inside);
+}

include/igl/cgal/is_inside.h

@@ -0,0 +1,82 @@
+#ifndef IS_INSIDE
+#define IS_INSIDE
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+    namespace cgal {
+
+        // Determine if mesh (V1, F1) is inside of mesh (V2, F2).
+        //
+        // Precondition:
+        // Both mesh must represent closed, self-intersection free,
+        // non-degenerated surfaces that
+        // are the boundary of 3D volumes. In addition, (V1, F1) must not
+        // intersect with (V2, F2).
+        //
+        // Inputs:
+        //   V1  #V1 by 3 list of vertex position of mesh 1
+        //   F1  #F1 by 3 list of triangles indices into V1
+        //   I1  #I1 list of indices into F1, faces to consider.
+        //   V2  #V2 by 3 list of vertex position of mesh 2
+        //   F2  #F2 by 3 list of triangles indices into V2
+        //   I2  #I2 list of indices into F2, faces to consider.
+        //
+        // Outputs:
+        //   return true iff (V1, F1) is entirely inside of (V2, F2).
+        template<typename DerivedV, typename DerivedF, typename DerivedI>
+            IGL_INLINE bool is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V1,
+                    const Eigen::PlainObjectBase<DerivedF>& F1,
+                    const Eigen::PlainObjectBase<DerivedI>& I1,
+                    const Eigen::PlainObjectBase<DerivedV>& V2,
+                    const Eigen::PlainObjectBase<DerivedF>& F2,
+                    const Eigen::PlainObjectBase<DerivedI>& I2);
+
+        template<typename DerivedV, typename DerivedF>
+            IGL_INLINE bool is_inside(
+                    const Eigen::PlainObjectBase<DerivedV>& V1,
+                    const Eigen::PlainObjectBase<DerivedF>& F1,
+                    const Eigen::PlainObjectBase<DerivedV>& V2,
+                    const Eigen::PlainObjectBase<DerivedF>& F2);
+
+        // Determine if queries points are inside of mesh (V, F).
+        //
+        // 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.
+        //
+        // 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 is_inside_batch(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedI>& I,
+                    const Eigen::PlainObjectBase<DerivedP>& P,
+                    Eigen::PlainObjectBase<DerivedB>& inside);
+
+        template<typename DerivedV, typename DerivedF, typename DerivedP,
+            typename DerivedB>
+            IGL_INLINE void is_inside_batch(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedP>& P,
+                    Eigen::PlainObjectBase<DerivedB>& inside);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "is_inside.cpp"
+#endif
+#endif

include/igl/cgal/order_facets_around_edge.cpp

@@ -3,353 +3,382 @@
 
 namespace igl
 {
-  namespace cgal
-  {
-    namespace order_facets_around_edges_helper
+    namespace cgal
     {
-      template<typename T>
-        std::vector<size_t> index_sort(const std::vector<T>& data)
+        namespace order_facets_around_edges_helper
         {
-          const size_t len = data.size();
-          std::vector<size_t> order(len);
-          for (size_t i=0; i<len; i++)
-          {
-            order[i] = i;
-          }
-          auto comp = [&](size_t i, size_t j)
-          {
-            return data[i] < data[j];
-          };
-          std::sort(order.begin(), order.end(), comp);
-          return order;
+            template<typename T>
+            std::vector<size_t> index_sort(const std::vector<T>& data)
+            {
+                const size_t len = data.size();
+                std::vector<size_t> order(len);
+                for (size_t i=0; i<len; i++) 
+                {
+                    order[i] = i;
+                }
+                auto comp = [&](size_t i, size_t j) 
+                {
+                    return data[i] < data[j];
+                };
+                std::sort(order.begin(), order.end(), comp);
+                return order;
+            }
         }
     }
-  }
 }
 
 // adj_faces contains signed index starting from +- 1.
 template<
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedI >
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI >
 void igl::cgal::order_facets_around_edge(
-  const Eigen::PlainObjectBase<DerivedV>& V,
-  const Eigen::PlainObjectBase<DerivedF>& F,
-  size_t s, 
-  size_t d, 
-  const std::vector<int>& adj_faces,
-  Eigen::PlainObjectBase<DerivedI>& order)
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    size_t s,
+    size_t d, 
+    const std::vector<int>& adj_faces,
+    Eigen::PlainObjectBase<DerivedI>& order)
 {
-  using namespace igl::cgal::order_facets_around_edges_helper;
-
-  // Although we only need exact predicates in the algorithm,
-  // exact constructions are needed to avoid degeneracies due to
-  // casting to double.
-  typedef CGAL::Exact_predicates_exact_constructions_kernel K;
-  typedef K::Point_3 Point_3;
-  typedef K::Plane_3 Plane_3;
-
-  auto get_face_index = [&](int adj_f)->size_t
-  {
-    return abs(adj_f) - 1;
-  };
-
-  auto get_opposite_vertex = [&](size_t fid)->size_t
-  {
-    typedef typename DerivedF::Scalar Index;
-    if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
-    if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
-    if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
-    assert(false);
-    return -1;
-  };
-
-  // Handle base cases
-  if (adj_faces.size() == 0) 
-  {
-    order.resize(0, 1);
-    return;
-  } else if (adj_faces.size() == 1)
-  {
-    order.resize(1, 1);
-    order(0, 0) = 0;
-    return;
-  } else if (adj_faces.size() == 2)
-  {
-    const size_t o1 = get_opposite_vertex(get_face_index(adj_faces[0]));
-    const size_t o2 = get_opposite_vertex(get_face_index(adj_faces[1]));
-    const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
-    const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
-    const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
-    const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
-    order.resize(2, 1);
-    switch (CGAL::orientation(ps, pd, p1, p2))
+    using namespace igl::cgal::order_facets_around_edges_helper;
+
+    // Although we only need exact predicates in the algorithm,
+    // exact constructions are needed to avoid degeneracies due to
+    // casting to double.
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef K::Point_3 Point_3;
+    typedef K::Plane_3 Plane_3;
+
+    auto get_face_index = [&](int adj_f)->size_t
+    {
+        return abs(adj_f) - 1;
+    };
+
+    auto get_opposite_vertex = [&](size_t fid)->size_t
+    {
+        typedef typename DerivedF::Scalar Index;
+        if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
+        if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
+        if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
+        assert(false);
+        return -1;
+    };
+
+    // Handle base cases
+    if (adj_faces.size() == 0) 
+    {
+        order.resize(0, 1);
+        return;
+    } else if (adj_faces.size() == 1)
     {
-      case CGAL::POSITIVE:
-        order(0, 0) = 1;
-        order(1, 0) = 0;
-        break;
-      case CGAL::NEGATIVE:
+        order.resize(1, 1);
         order(0, 0) = 0;
-        order(1, 0) = 1;
-        break;
-      case CGAL::COPLANAR:
+        return;
+    } else if (adj_faces.size() == 2)
+    {
+        const size_t o1 = get_opposite_vertex(get_face_index(adj_faces[0]));
+        const size_t o2 = get_opposite_vertex(get_face_index(adj_faces[1]));
+        const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
+        const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
+        const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
+        const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
+        order.resize(2, 1);
+        switch (CGAL::orientation(ps, pd, p1, p2))
         {
-          Plane_3 P1(ps, pd, p1);
-          Plane_3 P2(ps, pd, p2);
-          if (P1.orthogonal_direction() == P2.orthogonal_direction()){
-            // Duplicated face, use index to break tie.
-            order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
-            order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
-          } else {
-            // Coplanar faces, one on each side of the edge.
-            // It is equally valid to order them (0, 1) or (1, 0).
-            // I cannot think of any reason to prefer one to the
-            // other.  So just use (0, 1) ordering by default.
-            order(0, 0) = 0;
-            order(1, 0) = 1;
-          }
+            case CGAL::POSITIVE:
+                order(0, 0) = 1;
+                order(1, 0) = 0;
+                break;
+            case CGAL::NEGATIVE:
+                order(0, 0) = 0;
+                order(1, 0) = 1;
+                break;
+            case CGAL::COPLANAR:
+                {
+                    switch (CGAL::coplanar_orientation(ps, pd, p1, p2)) {
+                        case CGAL::POSITIVE:
+                            // Duplicated face, use index to break tie.
+                            order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
+                            order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+                            break;
+                        case CGAL::NEGATIVE:
+                            // Coplanar faces, one on each side of the edge.
+                            // It is equally valid to order them (0, 1) or (1, 0).
+                            // I cannot think of any reason to prefer one to the
+                            // other.  So just use (0, 1) ordering by default.
+                            order(0, 0) = 0;
+                            order(1, 0) = 1;
+                            break;
+                        case CGAL::COLLINEAR:
+                            std::cerr << "Degenerated triangle detected." <<
+                                std::endl;
+                            assert(false);
+                            break;
+                        default:
+                            assert(false);
+                    }
+                }
+                break;
+            default:
+                assert(false);
         }
-        break;
-      default:
-        assert(false);
+        return;
     }
-    return;
-  }
-
-  const size_t num_adj_faces = adj_faces.size();
-  const size_t o = get_opposite_vertex( get_face_index(adj_faces[0]));
-  const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
-  const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
-  const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
-  const Plane_3 separator(p_s, p_d, p_o);
-  assert(!separator.is_degenerate());
-
-  std::vector<Point_3> opposite_vertices;
-  for (size_t i=0; i<num_adj_faces; i++) 
-  {
-    const size_t o = get_opposite_vertex( get_face_index(adj_faces[i]));
-    opposite_vertices.emplace_back(
-        V(o, 0), V(o, 1), V(o, 2));
-  }
-
-  std::vector<int> positive_side;
-  std::vector<int> negative_side;
-  std::vector<int> tie_positive_oriented;
-  std::vector<int> tie_negative_oriented;
-
-  std::vector<size_t> positive_side_index;
-  std::vector<size_t> negative_side_index;
-  std::vector<size_t> tie_positive_oriented_index;
-  std::vector<size_t> tie_negative_oriented_index;
-
-  for (size_t i=0; i<num_adj_faces; i++)
-  {
-    const int f = adj_faces[i];
-    const Point_3& p_a = opposite_vertices[i];
-    auto orientation = separator.oriented_side(p_a);
-    switch (orientation) {
-      case CGAL::ON_POSITIVE_SIDE:
-        positive_side.push_back(f);
-        positive_side_index.push_back(i);
-        break;
-      case CGAL::ON_NEGATIVE_SIDE:
-        negative_side.push_back(f);
-        negative_side_index.push_back(i);
-        break;
-      case CGAL::ON_ORIENTED_BOUNDARY:
-        {
-          const Plane_3 other(p_s, p_d, p_a);
-          const auto target_dir = separator.orthogonal_direction();
-          const auto query_dir = other.orthogonal_direction();
-          if (target_dir == query_dir) {
-            tie_positive_oriented.push_back(f);
-            tie_positive_oriented_index.push_back(i);
-          } else if (target_dir == -query_dir) {
-            tie_negative_oriented.push_back(f);
-            tie_negative_oriented_index.push_back(i);
-          } else {
-            assert(false);
-          }
+
+    const size_t num_adj_faces = adj_faces.size();
+    const size_t o = get_opposite_vertex( get_face_index(adj_faces[0]));
+    const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
+    const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
+    const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
+    const Plane_3 separator(p_s, p_d, p_o);
+    assert(!separator.is_degenerate());
+
+    std::vector<Point_3> opposite_vertices;
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
+        const size_t o = get_opposite_vertex( get_face_index(adj_faces[i]));
+        opposite_vertices.emplace_back(
+                V(o, 0), V(o, 1), V(o, 2));
+    }
+
+    std::vector<int> positive_side;
+    std::vector<int> negative_side;
+    std::vector<int> tie_positive_oriented;
+    std::vector<int> tie_negative_oriented;
+
+    std::vector<size_t> positive_side_index;
+    std::vector<size_t> negative_side_index;
+    std::vector<size_t> tie_positive_oriented_index;
+    std::vector<size_t> tie_negative_oriented_index;
+
+    for (size_t i=0; i<num_adj_faces; i++)
+    {
+        const int f = adj_faces[i];
+        const Point_3& p_a = opposite_vertices[i];
+        auto orientation = separator.oriented_side(p_a);
+        switch (orientation) {
+            case CGAL::ON_POSITIVE_SIDE:
+                positive_side.push_back(f);
+                positive_side_index.push_back(i);
+                break;
+            case CGAL::ON_NEGATIVE_SIDE:
+                negative_side.push_back(f);
+                negative_side_index.push_back(i);
+                break;
+            case CGAL::ON_ORIENTED_BOUNDARY:
+                {
+                    auto inplane_orientation = CGAL::coplanar_orientation(
+                            p_s, p_d, p_o, p_a);
+                    switch (inplane_orientation) {
+                        case CGAL::POSITIVE:
+                            tie_positive_oriented.push_back(f);
+                            tie_positive_oriented_index.push_back(i);
+                            break;
+                        case CGAL::NEGATIVE:
+                            tie_negative_oriented.push_back(f);
+                            tie_negative_oriented_index.push_back(i);
+                            break;
+                        case CGAL::COLLINEAR:
+                        default:
+                            assert(false);
+                            break;
+                    }
+                }
+                break;
+            default:
+                // Should not be here.
+                assert(false);
         }
-        break;
-      default:
-        // Should not be here.
-        assert(false);
     }
-  }
-
-  Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
-  order_facets_around_edge(V, F, s, d, positive_side, positive_order);
-  order_facets_around_edge(V, F, s, d, negative_side, negative_order);
-  std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
-  std::vector<size_t> tie_negative_order = index_sort(tie_negative_oriented);
-
-  // Copy results into order vector.
-  const size_t tie_positive_size = tie_positive_oriented.size();
-  const size_t tie_negative_size = tie_negative_oriented.size();
-  const size_t positive_size = positive_order.size();
-  const size_t negative_size = negative_order.size();
-
-  order.resize(
-    tie_positive_size + positive_size + tie_negative_size + negative_size,1);
-
-  size_t count=0;
-  for (size_t i=0; i<tie_positive_size; i++)
-  {
-    order(count+i, 0) = tie_positive_oriented_index[tie_positive_order[i]];
-  }
-  count += tie_positive_size;
-
-  for (size_t i=0; i<negative_size; i++)
-  {
-    order(count+i, 0) = negative_side_index[negative_order(i, 0)];
-  }
-  count += negative_size;
-
-  for (size_t i=0; i<tie_negative_size; i++)
-  {
-    order(count+i, 0) = tie_negative_oriented_index[tie_negative_order[i]];
-  }
-  count += tie_negative_size;
-
-  for (size_t i=0; i<positive_size; i++)
-  {
-    order(count+i, 0) = positive_side_index[positive_order(i, 0)];
-  }
-  count += positive_size;
-  assert(count == num_adj_faces);
-
-  // Find the correct start point.
-  size_t start_idx = 0;
-  for (size_t i=0; i<num_adj_faces; i++)
-  {
-    const Point_3& p_a = opposite_vertices[order(i, 0)];
-    const Point_3& p_b =
-      opposite_vertices[order((i+1)%num_adj_faces, 0)];
-    auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
-    if (orientation == CGAL::POSITIVE)
+
+    Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
+    order_facets_around_edge(V, F, s, d, positive_side, positive_order);
+    order_facets_around_edge(V, F, s, d, negative_side, negative_order);
+    std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
+    std::vector<size_t> tie_negative_order = index_sort(tie_negative_oriented);
+
+    // Copy results into order vector.
+    const size_t tie_positive_size = tie_positive_oriented.size();
+    const size_t tie_negative_size = tie_negative_oriented.size();
+    const size_t positive_size = positive_order.size();
+    const size_t negative_size = negative_order.size();
+
+    order.resize(
+      tie_positive_size + positive_size + tie_negative_size + negative_size,1);
+
+    size_t count=0;
+    for (size_t i=0; i<tie_positive_size; i++)
+    {
+        order(count+i, 0) = tie_positive_oriented_index[tie_positive_order[i]];
+    }
+    count += tie_positive_size;
+
+    for (size_t i=0; i<negative_size; i++) 
+    {
+        order(count+i, 0) = negative_side_index[negative_order(i, 0)];
+    }
+    count += negative_size;
+
+    for (size_t i=0; i<tie_negative_size; i++)
+    {
+        order(count+i, 0) = tie_negative_oriented_index[tie_negative_order[i]];
+    }
+    count += tie_negative_size;
+
+    for (size_t i=0; i<positive_size; i++)
+    {
+        order(count+i, 0) = positive_side_index[positive_order(i, 0)];
+    }
+    count += positive_size;
+    assert(count == num_adj_faces);
+
+    // Find the correct start point.
+    size_t start_idx = 0;
+    for (size_t i=0; i<num_adj_faces; i++)
     {
-      // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
-      // more than 180 degrees.
-      start_idx = (i+1)%num_adj_faces;
-      break;
-    } else if (orientation == CGAL::COPLANAR &&
-        Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
-        Plane_3(p_s, p_d, p_b).orthogonal_direction())
+        const Point_3& p_a = opposite_vertices[order(i, 0)];
+        const Point_3& p_b =
+            opposite_vertices[order((i+1)%num_adj_faces, 0)];
+        auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
+        if (orientation == CGAL::POSITIVE)
+        {
+            // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
+            // more than 180 degrees.
+            start_idx = (i+1)%num_adj_faces;
+            break;
+        } else if (orientation == CGAL::COPLANAR &&
+                Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
+                Plane_3(p_s, p_d, p_b).orthogonal_direction())
+        {
+            // All 4 points are coplanar, but p_a and p_b are on each side of
+            // the edge (p_s, p_d).  This means the angle between triangle
+            // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
+            start_idx = (i+1)%num_adj_faces;
+            break;
+        }
+    }
+    DerivedI circular_order = order;
+    for (size_t i=0; i<num_adj_faces; i++)
     {
-      // All 4 points are coplanar, but p_a and p_b are on each side of
-      // the edge (p_s, p_d).  This means the angle between triangle
-      // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
-      start_idx = (i+1)%num_adj_faces;
-      break;
+        order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
     }
-  }
-  DerivedI circular_order = order;
-  for (size_t i=0; i<num_adj_faces; i++)
-  {
-    order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
-  }
 }
 
 template<
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedI>
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI>
 IGL_INLINE
 void igl::cgal::order_facets_around_edge(
-  const Eigen::PlainObjectBase<DerivedV>& V,
-  const Eigen::PlainObjectBase<DerivedF>& F,
-  size_t s, 
-  size_t d, 
-  const std::vector<int>& adj_faces,
-  const Eigen::PlainObjectBase<DerivedV>& pivot_point,
-  Eigen::PlainObjectBase<DerivedI>& order) 
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        size_t s,
+        size_t d, 
+        const std::vector<int>& adj_faces,
+        const Eigen::PlainObjectBase<DerivedV>& pivot_point,
+        Eigen::PlainObjectBase<DerivedI>& order)
 {
 
-  assert(V.cols() == 3);
-  assert(F.cols() == 3);
-  auto signed_index_to_index = [&](int signed_idx)
-  {
-    return abs(signed_idx) -1;
-  };
-  auto get_opposite_vertex_index = [&](size_t fid)
-  {
-    typedef typename DerivedF::Scalar Index;
-    if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
-    if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
-    if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
-    assert(false);
-    // avoid warning
-    return -1;
-  };
-
-  const typename DerivedF::Scalar N = adj_faces.size();
-  const size_t num_faces = N + 1; // N adj faces + 1 pivot face
-
-  DerivedV vertices(num_faces + 2, 3);
-  for (size_t i=0; i<(size_t)N; i++)
-  {
-    const size_t fid = signed_index_to_index(adj_faces[i]);
-    vertices.row(i) = V.row(get_opposite_vertex_index(fid));
-  }
-  vertices.row(N  ) = pivot_point;
-  vertices.row(N+1) = V.row(s);
-  vertices.row(N+2) = V.row(d);
-
-  DerivedF faces(num_faces, 3);
-  for (size_t i=0; i<(size_t)N; i++) 
-  {
-    if (adj_faces[i] < 0)
+    assert(V.cols() == 3);
+    assert(F.cols() == 3);
+    auto signed_index_to_index = [&](int signed_idx)
+    {
+        return abs(signed_idx) -1;
+    };
+    auto get_opposite_vertex_index = [&](size_t fid)
+    {
+        typedef typename DerivedF::Scalar Index;
+        if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
+        if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
+        if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
+        assert(false);
+        // avoid warning
+        return -1;
+    };
+
     {
-      faces(i,0) = N+1; // s
-      faces(i,1) = N+2; // d
-      faces(i,2) = i  ;
-    } else {
-      faces(i,0) = N+2; // d
-      faces(i,1) = N+1; // s
-      faces(i,2) = i  ;
+        // Check if s, d and pivot are collinear.
+        typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+        K::Point_3 ps(V(s,0), V(s,1), V(s,2));
+        K::Point_3 pd(V(d,0), V(d,1), V(d,2));
+        K::Point_3 pp(pivot_point(0,0), pivot_point(0,1), pivot_point(0,2));
+        if (CGAL::collinear(ps, pd, pp)) {
+            throw "Pivot point is collinear with the outer edge!";
+        }
     }
-  }
-  // Last face is the pivot face.
-  faces(N, 0) = N+1;
-  faces(N, 1) = N+2;
-  faces(N, 2) = N;
-
-  std::vector<int> adj_faces_with_pivot(num_faces);
-  for (size_t i=0; i<num_faces; i++)
-  {
-    if (faces(i,0) == N+1 && faces(i,1) == N+2)
+
+    const size_t N = adj_faces.size();
+    const size_t num_faces = N + 1; // N adj faces + 1 pivot face
+
+    // Because face indices are used for tie breaking, the original face indices
+    // in the new faces array must be ascending.
+    auto comp = [&](int i, int j) {
+        return signed_index_to_index(i) < signed_index_to_index(j);
+    };
+    std::vector<int> ordered_adj_faces(adj_faces);
+    std::sort(ordered_adj_faces.begin(), ordered_adj_faces.end(), comp);
+
+    DerivedV vertices(num_faces + 2, 3);
+    for (size_t i=0; i<N; i++) {
+        const size_t fid = signed_index_to_index(ordered_adj_faces[i]);
+        vertices.row(i) = V.row(get_opposite_vertex_index(fid));
+    }
+    vertices.row(N  ) = pivot_point;
+    vertices.row(N+1) = V.row(s);
+    vertices.row(N+2) = V.row(d);
+
+    DerivedF faces(num_faces, 3);
+    for (size_t i=0; i<N; i++)
     {
-      adj_faces_with_pivot[i] = int(i+1) * -1;
-    } else
+        if (ordered_adj_faces[i] < 0) {
+            faces(i,0) = N+1; // s
+            faces(i,1) = N+2; // d
+            faces(i,2) = i  ;
+        } else {
+            faces(i,0) = N+2; // d
+            faces(i,1) = N+1; // s
+            faces(i,2) = i  ;
+        }
+    }
+    // Last face is the pivot face.
+    faces(N, 0) = N+1;
+    faces(N, 1) = N+2;
+    faces(N, 2) = N;
+
+    std::vector<int> adj_faces_with_pivot(num_faces);
+    for (size_t i=0; i<num_faces; i++)
     {
-      adj_faces_with_pivot[i] = int(i+1);
+        if (faces(i,0) == N+1 && faces(i,1) == N+2)
+        {
+            adj_faces_with_pivot[i] = int(i+1) * -1;
+        } else
+        {
+            adj_faces_with_pivot[i] = int(i+1);
+        }
     }
-  }
-
-  DerivedI order_with_pivot;
-  igl::cgal::order_facets_around_edge(
-      vertices, faces, N+1, N+2,
-      adj_faces_with_pivot, order_with_pivot);
-
-  assert(order_with_pivot.size() == num_faces);
-  order.resize(N);
-  size_t pivot_index = num_faces + 1;
-  for (size_t i=0; i<num_faces; i++)
-  {
-    if (order_with_pivot[i] == N)
+
+    DerivedI order_with_pivot;
+    igl::cgal::order_facets_around_edge(
+            vertices, faces, N+1, N+2,
+            adj_faces_with_pivot, order_with_pivot);
+
+    assert(order_with_pivot.size() == num_faces);
+    order.resize(N);
+    size_t pivot_index = num_faces + 1;
+    for (size_t i=0; i<num_faces; i++)
     {
-      pivot_index = i;
-      break;
+        if (order_with_pivot[i] == N)
+        {
+            pivot_index = i;
+            break;
+        }
     }
-  }
-  assert(pivot_index < num_faces);
+    assert(pivot_index < num_faces);
 
-  for (size_t i=0; i<(size_t)N; i++)
-  {
-    order[i] = order_with_pivot[(pivot_index+i+1)%num_faces];
-  }
+    for (size_t i=0; i<N; i++)
+    {
+        order[i] = order_with_pivot[(pivot_index+i+1)%num_faces];
+    }
 }
 
 

include/igl/cgal/outer_hull.cpp

@@ -5,6 +5,7 @@
 // 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 "is_inside.h"
 #include "outer_hull.h"
 #include "order_facets_around_edges.h"
 #include "outer_facet.h"
@@ -17,6 +18,7 @@
 #include "../per_face_normals.h"
 #include "../writePLY.h"
 #include "../sort_angles.h"
+#include "../writePLY.h"
 
 #include <Eigen/Geometry>
 #include <vector>
@@ -325,19 +327,19 @@ IGL_INLINE void igl::cgal::outer_hull(
 
   // Is A inside B? Assuming A and B are consistently oriented but closed and
   // non-intersecting.
-  const auto & is_component_inside_other = [](
-    const Eigen::MatrixXd & V,
+  const auto & has_overlapping_bbox = [](
+    const Eigen::PlainObjectBase<DerivedV> & V,
     const MatrixXV & BC,
     const MatrixXG & A,
     const MatrixXJ & AJ,
     const MatrixXG & B)->bool
   {
     const auto & bounding_box = [](
-      const Eigen::MatrixXd & V,
+      const Eigen::PlainObjectBase<DerivedV> & V,
       const MatrixXG & F)->
-      Eigen::MatrixXd
+        DerivedV
     {
-      Eigen::MatrixXd BB(2,3);
+      DerivedV BB(2,3);
       BB<<
          1e26,1e26,1e26,
         -1e26,-1e26,-1e26;
@@ -346,75 +348,73 @@ IGL_INLINE void igl::cgal::outer_hull(
       {
         for(size_t c = 0;c<3;c++)
         {
-          const auto & vfc = V.row(F(f,c));
-          BB.row(0) = BB.row(0).array().min(vfc.array()).eval();
-          BB.row(1) = BB.row(1).array().max(vfc.array()).eval();
+          const auto & vfc = V.row(F(f,c)).eval();
+          BB(0,0) = std::min(BB(0,0), vfc(0,0));
+          BB(0,1) = std::min(BB(0,1), vfc(0,1));
+          BB(0,2) = std::min(BB(0,2), vfc(0,2));
+          BB(1,0) = std::max(BB(1,0), vfc(0,0));
+          BB(1,1) = std::max(BB(1,1), vfc(0,1));
+          BB(1,2) = std::max(BB(1,2), vfc(0,2));
         }
       }
       return BB;
     };
     // A lot of the time we're dealing with unrelated, distant components: cull
     // them.
-    Eigen::MatrixXd ABB = bounding_box(V,A);
-    Eigen::MatrixXd BBB = bounding_box(V,B);
+    DerivedV ABB = bounding_box(V,A);
+    DerivedV BBB = bounding_box(V,B);
     if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
         (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
     {
       // bounding boxes do not overlap
       return false;
+    } else {
+      return true;
     }
-    ////////////////////////////////////////////////////////////////////////
-    // POTENTIAL ROBUSTNESS WEAK AREA
-    ////////////////////////////////////////////////////////////////////////
-    //
-
-    // winding_number_3 expects colmajor
-    // q could be so close (<~1e-15) to B that the winding number is not a robust way to
-    // determine inside/outsideness. We could try to find a _better_ q which is
-    // farther away, but couldn't they all be bad?
-    double q[3] = {
-        CGAL::to_double(BC(AJ(0), 0)),
-        CGAL::to_double(BC(AJ(0), 1)),
-        CGAL::to_double(BC(AJ(0), 2)) };
-    // In a perfect world, it's enough to test a single point.
-    double w;
-    winding_number_3(
-      V.data(),V.rows(),
-      B.data(),B.rows(),
-      q,1,&w);
-    return w > 0.5 || w < -0.5;
   };
 
-  Eigen::MatrixXd Vcol(V.rows(), V.cols());
-  for (size_t i=0; i<(size_t)V.rows(); i++) {
-      for (size_t j=0; j<(size_t)V.cols(); j++) {
-          Vcol(i, j) = CGAL::to_double(V(i, j));
-      }
-  }
-
   // Reject components which are completely inside other components
   vector<bool> keep(ncc,true);
   size_t nG = 0;
   // This is O( ncc * ncc * m)
   for(size_t id = 0;id<ncc;id++)
   {
+    if (!keep[id]) continue;
+    std::vector<size_t> unresolved;
     for(size_t oid = 0;oid<ncc;oid++)
     {
-      if(id == oid)
+      if(id == oid || !keep[oid])
       {
         continue;
       }
-      const bool inside = is_component_inside_other(Vcol,BC,vG[id],vJ[id],vG[oid]);
-#ifdef IGL_OUTER_HULL_DEBUG
-      cout<<id<<" is inside "<<oid<<" ? "<<inside<<endl;
-#endif
-      keep[id] = keep[id] && !inside;
+      if (has_overlapping_bbox(V, BC, vG[id], vJ[id], vG[oid])) {
+          unresolved.push_back(oid);
+      }
     }
-    if(keep[id])
-    {
-      nG += vJ[id].rows();
+    const size_t num_unresolved_components = unresolved.size();
+    DerivedV query_points(num_unresolved_components, 3);
+    for (size_t i=0; i<num_unresolved_components; i++) {
+        const size_t oid = unresolved[i];
+        DerivedF f = vG[oid].row(0);
+        query_points(i,0) = (V(f(0,0), 0) + V(f(0,1), 0) + V(f(0,2), 0))/3.0;
+        query_points(i,1) = (V(f(0,0), 1) + V(f(0,1), 1) + V(f(0,2), 1))/3.0;
+        query_points(i,2) = (V(f(0,0), 2) + V(f(0,1), 2) + V(f(0,2), 2))/3.0;
+    }
+    Eigen::VectorXi inside;
+    igl::cgal::is_inside_batch(V, vG[id], query_points, inside);
+    assert(inside.size() == num_unresolved_components);
+    for (size_t i=0; i<num_unresolved_components; i++) {
+        if (inside[i]) {
+            const size_t oid = unresolved[i];
+            keep[oid] = false;
+        }
     }
   }
+  for (size_t id = 0; id<ncc; id++) {
+      if (keep[id]) {
+          nG += vJ[id].rows();
+      }
+  }
 
   // collect G and J across components
   G.resize(nG,3);

include/igl/outer_element.cpp

@@ -23,59 +23,59 @@ IGL_INLINE void igl::outer_vertex(
         IndexType & v_index,
         Eigen::PlainObjectBase<DerivedA> & A)
 {
-  // Algorithm: 
-  //    Find an outer vertex (i.e. vertex reachable from infinity)
-  //    Return the vertex with the largest X value.
-  //    If there is a tie, pick the one with largest Y value.
-  //    If there is still a tie, pick the one with the largest Z value.
-  //    If there is still a tie, then there are duplicated vertices within the
-  //    mesh, which violates the precondition.
-  typedef typename DerivedF::Scalar Index;
-  const Index INVALID = std::numeric_limits<Index>::max();
-  const size_t num_selected_faces = I.rows();
-  std::vector<size_t> candidate_faces;
-  Index outer_vid = INVALID;
-  typename DerivedV::Scalar outer_val = 0;
-  for (size_t i=0; i<num_selected_faces; i++)
-  {
-    size_t f = I(i);
-    for (size_t j=0; j<3; j++)
+    // Algorithm: 
+    //    Find an outer vertex (i.e. vertex reachable from infinity)
+    //    Return the vertex with the largest X value.
+    //    If there is a tie, pick the one with largest Y value.
+    //    If there is still a tie, pick the one with the largest Z value.
+    //    If there is still a tie, then there are duplicated vertices within the
+    //    mesh, which violates the precondition.
+    typedef typename DerivedF::Scalar Index;
+    const Index INVALID = std::numeric_limits<Index>::max();
+    const size_t num_selected_faces = I.rows();
+    std::vector<size_t> candidate_faces;
+    Index outer_vid = INVALID;
+    typename DerivedV::Scalar outer_val = 0;
+    for (size_t i=0; i<num_selected_faces; i++)
     {
-      Index v = F(f, j);
-      auto vx = V(v, 0);
-      if (outer_vid == INVALID || vx > outer_val)
-      {
-        outer_val = vx;
-        outer_vid = v;
-        candidate_faces = {f};
-      } else if (v == outer_vid)
-      {
-        candidate_faces.push_back(f);
-      } else if (vx == outer_val)
-      {
-        // Break tie.
-        auto vy = V(v,1);
-        auto vz = V(v, 2);
-        auto outer_y = V(outer_vid, 1);
-        auto outer_z = V(outer_vid, 2);
-        assert(!(vy == outer_y && vz == outer_z));
-        bool replace = (vy > outer_y) ||
-          ((vy == outer_y) && (vz > outer_z));
-        if (replace)
+        size_t f = I(i);
+        for (size_t j=0; j<3; j++)
         {
-          outer_val = vx;
-          outer_vid = v;
-          candidate_faces = {f};
+            Index v = F(f, j);
+            auto vx = V(v, 0);
+            if (outer_vid == INVALID || vx > outer_val)
+            {
+                outer_val = vx;
+                outer_vid = v;
+                candidate_faces = {f};
+            } else if (v == outer_vid)
+            {
+                candidate_faces.push_back(f);
+            } else if (vx == outer_val)
+            {
+                // Break tie.
+                auto vy = V(v,1);
+                auto vz = V(v, 2);
+                auto outer_y = V(outer_vid, 1);
+                auto outer_z = V(outer_vid, 2);
+                assert(!(vy == outer_y && vz == outer_z));
+                bool replace = (vy > outer_y) ||
+                    ((vy == outer_y) && (vz > outer_z));
+                if (replace)
+                {
+                    outer_val = vx;
+                    outer_vid = v;
+                    candidate_faces = {f};
+                }
+            }
         }
-      }
     }
-  }
 
-  assert(outer_vid != INVALID);
-  assert(candidate_faces.size() > 0);
-  v_index = outer_vid;
-  A.resize(candidate_faces.size());
-  std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
+    assert(outer_vid != INVALID);
+    assert(candidate_faces.size() > 0);
+    v_index = outer_vid;
+    A.resize(candidate_faces.size());
+    std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
 }
 
 template<
@@ -86,94 +86,122 @@ template<
     typename DerivedA
     >
 IGL_INLINE void igl::outer_edge(
-       const Eigen::PlainObjectBase<DerivedV> & V,
-       const Eigen::PlainObjectBase<DerivedF> & F,
-       const Eigen::PlainObjectBase<DerivedI> & I,
-       IndexType & v1,
-       IndexType & v2,
-       Eigen::PlainObjectBase<DerivedA> & A) {
-   // Algorithm:
-   //    Find an outer vertex first.
-   //    Find the incident edge with largest slope when projected onto XY
-   //      plane.
-   //    If there is still a tie, break it using the projected slope onto ZX
-   //      plane.
-   //    If there is still a tie, then there are multiple overlapping edges,
-   //      which violates the precondition.
-   typedef typename DerivedV::Scalar Scalar;
-   typedef typename DerivedV::Index Index;
-   typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
-   typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
-   const Index INVALID = std::numeric_limits<Index>::max();
-
-   Index outer_vid;
-   Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
-   outer_vertex(V, F, I, outer_vid, candidate_faces);
-   const ScalarArray3& outer_v = V.row(outer_vid);
-   assert(candidate_faces.size() > 0);
-
-   auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
-   {
-     if (f[0] == vid) return 0;
-     if (f[1] == vid) return 1;
-     if (f[2] == vid) return 2;
-     assert(false);
-     return -1;
-   };
-
-   Scalar outer_slope_YX = 0;
-   Scalar outer_slope_ZX = 0;
-   Index outer_opp_vid = INVALID;
-   bool infinite_slope_detected = false;
-   std::vector<Index> incident_faces;
-   auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
-     if (opp_vid == outer_opp_vid)
-     {
-       incident_faces.push_back(fid);
-       return;
-     }
-
-     const ScalarArray3 opp_v = V.row(opp_vid);
-     if (!infinite_slope_detected && outer_v[0] != opp_v[0]) 
-     {
-       // Finite slope
-       const ScalarArray3 diff = opp_v - outer_v;
-       const Scalar slope_YX = diff[1] / diff[0];
-       const Scalar slope_ZX = diff[2] / diff[0];
-       if (outer_opp_vid == INVALID ||
-               slope_YX > outer_slope_YX ||
-               (slope_YX == outer_slope_YX &&
-                slope_ZX > outer_slope_ZX)) {
-           outer_opp_vid = opp_vid;
-           outer_slope_YX = slope_YX;
-           outer_slope_ZX = slope_ZX;
-           incident_faces = {fid};
-       }
-     } else if (!infinite_slope_detected)
-     {
-       // Infinite slope
-       outer_opp_vid = opp_vid;
-       infinite_slope_detected = true;
-       incident_faces = {fid};
-     }
-   };
-
-   const size_t num_candidate_faces = candidate_faces.size();
-   for (size_t i=0; i<num_candidate_faces; i++) 
-   {
-     const Index fid = candidate_faces(i);
-     const IndexArray3& f = F.row(fid);
-     size_t id = get_vertex_index(f, outer_vid);
-     Index next_vid = f((id+1)%3);
-     Index prev_vid = f((id+2)%3);
-     check_and_update_outer_edge(next_vid, fid);
-     check_and_update_outer_edge(prev_vid, fid);
-   }
-
-   v1 = outer_vid;
-   v2 = outer_opp_vid;
-   A.resize(incident_faces.size());
-   std::copy(incident_faces.begin(), incident_faces.end(), A.data());
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & v1,
+        IndexType & v2,
+        Eigen::PlainObjectBase<DerivedA> & A) {
+    // Algorithm:
+    //    Find an outer vertex first.
+    //    Find the incident edge with largest abs slope when projected onto XY plane.
+    //    If there is a tie, check the signed slope and use the positive one.
+    //    If there is still a tie, break it using the projected slope onto ZX plane.
+    //    If there is still a tie, again check the signed slope and use the positive one.
+    //    If there is still a tie, then there are multiple overlapping edges,
+    //    which violates the precondition.
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
+    typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
+    const Index INVALID = std::numeric_limits<Index>::max();
+
+    Index outer_vid;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
+    outer_vertex(V, F, I, outer_vid, candidate_faces);
+    const ScalarArray3& outer_v = V.row(outer_vid);
+    assert(candidate_faces.size() > 0);
+
+    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
+    {
+        if (f[0] == vid) return 0;
+        if (f[1] == vid) return 1;
+        if (f[2] == vid) return 2;
+        assert(false);
+        return -1;
+    };
+
+    auto unsigned_value = [](Scalar v) -> Scalar {
+        if (v < 0) return v * -1;
+        else return v;
+    };
+
+    Scalar outer_slope_YX = 0;
+    Scalar outer_slope_ZX = 0;
+    Index outer_opp_vid = INVALID;
+    bool infinite_slope_detected = false;
+    std::vector<Index> incident_faces;
+    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
+        if (opp_vid == outer_opp_vid)
+        {
+            incident_faces.push_back(fid);
+            return;
+        }
+
+        const ScalarArray3 opp_v = V.row(opp_vid);
+        if (!infinite_slope_detected && outer_v[0] != opp_v[0])
+        {
+            // Finite slope
+            const ScalarArray3 diff = opp_v - outer_v;
+            const Scalar slope_YX = diff[1] / diff[0];
+            const Scalar slope_ZX = diff[2] / diff[0];
+            const Scalar u_slope_YX = unsigned_value(slope_YX);
+            const Scalar u_slope_ZX = unsigned_value(slope_ZX);
+            bool update = false;
+            if (outer_opp_vid == INVALID) {
+                update = true;
+            } else {
+                const Scalar u_outer_slope_YX = unsigned_value(outer_slope_YX);
+                if (u_slope_YX > u_outer_slope_YX) {
+                    update = true;
+                } else if (u_slope_YX == u_outer_slope_YX &&
+                        slope_YX > outer_slope_YX) {
+                    update = true;
+                } else if (slope_YX == outer_slope_YX) {
+                    const Scalar u_outer_slope_ZX =
+                        unsigned_value(outer_slope_ZX);
+                    if (u_slope_ZX > u_outer_slope_ZX) {
+                        update = true;
+                    } else if (u_slope_ZX == u_outer_slope_ZX &&
+                            slope_ZX > outer_slope_ZX) {
+                        update = true;
+                    } else if (slope_ZX == u_outer_slope_ZX) {
+                        assert(false);
+                    }
+                }
+            }
+
+            if (update) {
+                outer_opp_vid = opp_vid;
+                outer_slope_YX = slope_YX;
+                outer_slope_ZX = slope_ZX;
+                incident_faces = {fid};
+            }
+        } else if (!infinite_slope_detected)
+        {
+            // Infinite slope
+            outer_opp_vid = opp_vid;
+            infinite_slope_detected = true;
+            incident_faces = {fid};
+        }
+    };
+
+    const size_t num_candidate_faces = candidate_faces.size();
+    for (size_t i=0; i<num_candidate_faces; i++)
+    {
+        const Index fid = candidate_faces(i);
+        const IndexArray3& f = F.row(fid);
+        size_t id = get_vertex_index(f, outer_vid);
+        Index next_vid = f((id+1)%3);
+        Index prev_vid = f((id+2)%3);
+        check_and_update_outer_edge(next_vid, fid);
+        check_and_update_outer_edge(prev_vid, fid);
+    }
+
+    v1 = outer_vid;
+    v2 = outer_opp_vid;
+    A.resize(incident_faces.size());
+    std::copy(incident_faces.begin(), incident_faces.end(), A.data());
 }
 
 template<