Refactor outer_hull and is_inside method.

Former-commit-id: 9a9a5a057c030c0110b7cdd017d4a63acff17c4a

include/igl/cgal/is_inside.cpp

@@ -266,67 +266,20 @@ IGL_INLINE bool igl::cgal::is_inside(
         const Eigen::PlainObjectBase<DerivedV>& V2,
         const Eigen::PlainObjectBase<DerivedF>& F2,
         const Eigen::PlainObjectBase<DerivedI>& I2) {
-    using namespace igl::cgal::is_inside_helper;
     assert(F1.rows() > 0);
     assert(I1.rows() > 0);
     assert(F2.rows() > 0);
     assert(I2.rows() > 0);
 
-    //assert(!intersect_each_other(V1, F1, I1, V2, F2, I2));
-
-    const size_t num_faces = I2.rows();
-    std::vector<Triangle> triangles;
-    for (size_t i=0; i<num_faces; i++) {
-        const Eigen::Vector3i f = F2.row(I2(i, 0));
-        triangles.emplace_back(
-                Point_3(V2(f[0], 0), V2(f[0], 1), V2(f[0], 2)),
-                Point_3(V2(f[1], 0), V2(f[1], 1), V2(f[1], 2)),
-                Point_3(V2(f[2], 0), V2(f[2], 1), V2(f[2], 2)));
-        assert(!triangles.back().is_degenerate());
-    }
-    Tree tree(triangles.begin(), triangles.end());
-    tree.accelerate_distance_queries();
-
     const Eigen::Vector3i& f = F1.row(I1(0, 0));
-    const Point_3 query(
+    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);
-    // Computing the closest point to mesh2 is the only exact construction
-    // needed in the algorithm.
-    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(
-                V2, F2, I2, fid, closest_point, element_index)) {
-        case VERTEX:
-            {
-                //std::cout << "vertex case" << std::endl;
-                const size_t s = F2(I2(fid, 0), element_index);
-                return determine_point_vertex_orientation(
-                        V2, F2, I2, query, s);
-            }
-            break;
-        case EDGE:
-            {
-                //std::cout << "edge case" << std::endl;
-                const size_t s = F2(I2(fid, 0), (element_index+1)%3);
-                const size_t d = F2(I2(fid, 0), (element_index+2)%3);
-                return determine_point_edge_orientation(
-                        V2, F2, I2, query, s, d);
-            }
-            break;
-        case FACE:
-            //std::cout << "face case" << std::endl;
-            return determine_point_face_orientation(V2, F2, I2, query, fid);
-            break;
-        default:
-            assert(false);
-    }
-    assert(false);
-    return false;
+    Eigen::VectorXi inside;
+    igl::cgal::is_inside(V2, F2, I2, query, inside);
+    assert(inside.size() == 1);
+    return inside[0];
 }
 
 template<typename DerivedV, typename DerivedF>
@@ -340,3 +293,75 @@ IGL_INLINE bool igl::cgal::is_inside(
     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(
+        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(
+        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(V, F, I, P, inside);
+}

include/igl/cgal/is_inside.h

@@ -40,6 +40,39 @@ namespace igl {
                     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(
+                    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(
+                    const Eigen::PlainObjectBase<DerivedV>& V,
+                    const Eigen::PlainObjectBase<DerivedF>& F,
+                    const Eigen::PlainObjectBase<DerivedP>& P,
+                    Eigen::PlainObjectBase<DerivedB>& inside);
     }
 }
 

include/igl/cgal/outer_hull.cpp

@@ -327,15 +327,15 @@ 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 DerivedV & 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 DerivedV & V,
+      const Eigen::PlainObjectBase<DerivedV> & V,
       const MatrixXG & F)->
         DerivedV
     {
@@ -349,8 +349,12 @@ 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();
+          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(0,0), vfc(0,0));
+          BB(1,1) = std::max(BB(0,1), vfc(0,1));
+          BB(1,2) = std::max(BB(0,2), vfc(0,2));
         }
       }
       return BB;
@@ -365,7 +369,7 @@ IGL_INLINE void igl::cgal::outer_hull(
       // bounding boxes do not overlap
       return false;
     } else {
-        return is_inside(V, A, V, B);
+      return true;
     }
   };
 
@@ -375,22 +379,41 @@ IGL_INLINE void igl::cgal::outer_hull(
   // 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;
       }
-      bool inside = is_component_inside_other(V,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(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