fix asymptotic performance bug in extract cells

Former-commit-id: 584d3983c915d1a8cd5cae4612eef09afece63bd

include/igl/copyleft/cgal/closest_facet.cpp

@@ -12,13 +12,8 @@
 #include <stdexcept>
 #include <unordered_map>
 
-#include <CGAL/AABB_tree.h>
-#include <CGAL/AABB_traits.h>
-#include <CGAL/AABB_triangle_primitive.h>
-#include <CGAL/intersections.h>
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-
 #include "order_facets_around_edge.h"
+#include "submesh_aabb_tree.h"
 #include "../../vertex_triangle_adjacency.h"
 #include "../../writePLY.h"
 
@@ -39,7 +34,9 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     const std::vector<std::vector<uE2EType> >& uE2E,
     const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
     Eigen::PlainObjectBase<DerivedR>& R,
-    Eigen::PlainObjectBase<DerivedS>& S) {
+    Eigen::PlainObjectBase<DerivedS>& S)
+{
+
   typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
   typedef Kernel::Point_3 Point_3;
   typedef Kernel::Plane_3 Plane_3;
@@ -55,29 +52,64 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
         "Closest facet cannot be computed on empty mesh.");
   }
 
-  std::vector<std::vector<size_t> > VF;
-  std::vector<std::vector<size_t> > VFi;
+  std::vector<std::vector<size_t> > VF, VFi;
   igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
+  std::vector<bool> in_I;
+  std::vector<Triangle> triangles;
+  Tree tree;
+  submesh_aabb_tree(V,F,I,tree,triangles,in_I);
+
+  return closest_facet(
+    V,F,I,P,uE2E,EMAP,VF,VFi,tree,triangles,in_I,R,S);
+}
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename Kernel,
+  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,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    const std::vector<std::vector<size_t> > & VF,
+    const std::vector<std::vector<size_t> > & VFi,
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<
+        Kernel, 
+        CGAL::AABB_triangle_primitive<
+          Kernel, typename std::vector<
+            typename Kernel::Triangle_3 >::iterator > > > & tree,
+    const std::vector<typename Kernel::Triangle_3 > & triangles,
+    const std::vector<bool> & in_I,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S)
+{
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename Kernel::Plane_3 Plane_3;
+  typedef typename Kernel::Segment_3 Segment_3;
+  typedef typename Kernel::Triangle_3 Triangle;
+  typedef typename std::vector<Triangle>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
 
-  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");
-    }
+  if (F.rows() <= 0 || I.rows() <= 0) {
+    throw std::runtime_error(
+        "Closest facet cannot be computed on empty mesh.");
   }
-  Tree tree(triangles.begin(), triangles.end());
-  tree.accelerate_distance_queries();
 
-  auto on_the_positive_side = [&](size_t fid, const Point_3& p) -> bool{
+  auto on_the_positive_side = [&](size_t fid, const Point_3& p) -> bool
+  {
     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));
@@ -99,7 +131,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     return false;
   };
 
-  auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
+  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;
@@ -143,23 +176,28 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
       size_t query_idx,
       const size_t s, const size_t d,
       size_t preferred_facet,
-      bool& orientation) -> size_t {
+      bool& orientation) -> size_t 
+  {
     Point_3 query_point(
-        P(query_idx, 0),
-        P(query_idx, 1),
-        P(query_idx, 2));
+      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))) {
+        (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))) {
+        (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))) {
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 1))) 
+    {
       corner_idx = 0;
-    } else {
+    } else 
+    {
       std::cerr << "s: " << s << "\t d:" << d << std::endl;
       std::cerr << F.row(preferred_facet) << std::endl;
       throw std::runtime_error(
@@ -169,9 +207,11 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     auto ueid = EMAP(preferred_facet + corner_idx * F.rows());
     auto eids = uE2E[ueid];
     std::vector<size_t> intersected_face_indices;
-    for (auto eid : eids) {
+    for (auto eid : eids) 
+    {
       const size_t fid = eid % F.rows();
-      if (in_I[fid]) {
+      if (in_I[fid]) 
+      {
         intersected_face_indices.push_back(fid);
       }
     }
@@ -188,7 +228,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
         });
 
     assert(num_intersected_faces >= 1);
-    if (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.
@@ -200,8 +241,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     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);
+      intersected_face_signed_indices,
+      pivot, order);
 
     // Although first and last are equivalent, make the choice based on
     // preferred_facet.
@@ -381,7 +422,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     bool fid_ori = false;
 
     // Gether all facets sharing the closest point.
-    std::vector<Tree::Primitive_id> intersected_faces;
+    typename std::vector<typename 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();
@@ -389,7 +430,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     std::transform(intersected_faces.begin(),
         intersected_faces.end(),
         intersected_face_indices.begin(),
-        [&](const Tree::Primitive_id& itr) -> size_t
+        [&](const typename Tree::Primitive_id& itr) -> size_t
         { return I(itr-triangles.begin(), 0); });
 
     size_t element_index;

include/igl/copyleft/cgal/closest_facet.h

@@ -13,6 +13,12 @@
 #include <Eigen/Core>
 #include <vector>
 
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 namespace igl
 {
   namespace copyleft
@@ -32,14 +38,14 @@ namespace igl
       //   S  #P list of bools indicating on which side of the closest facet
       //      each query point lies.
       template<
-          typename DerivedV,
-          typename DerivedF,
-          typename DerivedI,
-          typename DerivedP,
-          typename uE2EType,
-          typename DerivedEMAP,
-          typename DerivedR,
-          typename DerivedS >
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        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,
@@ -49,23 +55,51 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
-
       template<
-          typename DerivedV,
-          typename DerivedF,
-          typename DerivedP,
-          typename uE2EType,
-          typename DerivedEMAP,
-          typename DerivedR,
-          typename DerivedS >
+        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);
+        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);
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename DerivedP,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename Kernel,
+        typename DerivedR,
+        typename DerivedS >
+      IGL_INLINE void closest_facet(
+          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,
+          const std::vector<std::vector<size_t> > & VF,
+          const std::vector<std::vector<size_t> > & VFi,
+          const CGAL::AABB_tree<
+            CGAL::AABB_traits<
+              Kernel, 
+              CGAL::AABB_triangle_primitive<
+                Kernel, typename std::vector<
+                  typename Kernel::Triangle_3 >::iterator > > > & tree,
+          const std::vector<typename Kernel::Triangle_3 > & triangles,
+          const std::vector<bool> & in_I,
+          Eigen::PlainObjectBase<DerivedR>& R,
+          Eigen::PlainObjectBase<DerivedS>& S);
     }
   }
 }

include/igl/copyleft/cgal/extract_cells.cpp

@@ -7,14 +7,22 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 //
 #include "extract_cells.h"
+#include "closest_facet.h"
+#include "order_facets_around_edge.h"
+#include "outer_facet.h"
+#include "submesh_aabb_tree.h"
 #include "../../extract_manifold_patches.h"
 #include "../../facet_components.h"
+#include "../../get_seconds.h"
 #include "../../triangle_triangle_adjacency.h"
 #include "../../unique_edge_map.h"
-#include "../../get_seconds.h"
-#include "closest_facet.h"
-#include "order_facets_around_edge.h"
-#include "outer_facet.h"
+#include "../../vertex_triangle_adjacency.h"
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 #include <iostream>
 #include <vector>
@@ -75,6 +83,17 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
   const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
   Eigen::PlainObjectBase<DerivedC>& cells) 
 {
+
+  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;
+
 #ifdef EXTRACT_CELLS_DEBUG
   const auto & tictoc = []() -> double
   {
@@ -121,11 +140,25 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     components[C[i]].push_back(i);
   }
   // Convert vector lists to Eigen lists...
+  // and precompute data-structures for each component
+  std::vector<std::vector<size_t> > VF,VFi;
+  igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
   std::vector<Eigen::VectorXi> Is(num_components);
+  std::vector<
+    CGAL::AABB_tree<
+      CGAL::AABB_traits<
+        Kernel, 
+        CGAL::AABB_triangle_primitive<
+          Kernel, std::vector<
+            Kernel::Triangle_3 >::iterator > > > > trees(num_components);
+  std::vector< std::vector<Kernel::Triangle_3 > > 
+    triangle_lists(num_components);
+  std::vector<std::vector<bool> > in_Is(num_components);
   for (size_t i=0; i<num_components; i++)
   {
     Is[i].resize(components[i].size());
     std::copy(components[i].begin(), components[i].end(),Is[i].data());
+    submesh_aabb_tree(V,F,Is[i],trees[i],triangle_lists[i],in_Is[i]);
   }
 
   // Find outer facets, their orientations and cells for each component
@@ -219,9 +252,25 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
       // Gather closest facets in ith component to each query point and their
       // orientations
       const auto& I = Is[i];
+      const auto& tree = trees[i];
+      const auto& in_I = in_Is[i];
+      const auto& triangles = triangle_lists[i];
+
       Eigen::VectorXi closest_facets, closest_facet_orientations;
-      closest_facet(V, F, I, queries,
-        uE2E, EMAP, closest_facets, closest_facet_orientations);
+      closest_facet(
+        V,
+        F, 
+        I, 
+        queries,
+        uE2E, 
+        EMAP, 
+        VF,
+        VFi,
+        tree,
+        triangles,
+        in_I,
+        closest_facets, 
+        closest_facet_orientations);
       // Loop over all candidates
       for (size_t j=0; j<num_candidate_comps; j++)
       {

include/igl/copyleft/cgal/outer_hull.cpp

@@ -10,6 +10,12 @@
 #include "remesh_self_intersections.h"
 #include "../../remove_unreferenced.h"
 
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 template <
   typename DerivedV,
   typename DerivedF,
@@ -525,6 +531,4 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull_legacy(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::copyleft::cgal::outer_hull<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-template void igl::copyleft::cgal::outer_hull<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -74,7 +74,7 @@ IGL_INLINE size_t igl::copyleft::cgal::peel_outer_hull_layers(
       writePLY(ss.str(), vertices, Fr);
   }
 #endif
-    outer_hull(V,Fr,Fo,Jo,flipr);
+    outer_hull_legacy(V,Fr,Fo,Jo,flipr);
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
   writePLY(STR("outer-hull-output-"<<iter<<".ply"),V,Fo);
   cout<<"reindex, flip..."<<endl;

include/igl/copyleft/cgal/submesh_aabb_tree.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// 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 "submesh_aabb_tree.h"
+#include <stdexcept>
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename Kernel>
+IGL_INLINE void igl::copyleft::cgal::submesh_aabb_tree(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedI>& I,
+  CGAL::AABB_tree<
+    CGAL::AABB_traits<
+      Kernel, 
+      CGAL::AABB_triangle_primitive<
+        Kernel, typename std::vector<
+          typename Kernel::Triangle_3 >::iterator > > > & tree,
+  std::vector<typename Kernel::Triangle_3 > & triangles,
+  std::vector<bool> & in_I)
+{
+  in_I.resize(F.rows(), false);
+  const size_t num_faces = I.rows();
+  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(
+      typename Kernel::Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+      typename Kernel::Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+      typename Kernel::Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+#ifndef NDEBUG
+    if (triangles.back().is_degenerate()) 
+    {
+      throw std::runtime_error(
+        "Input facet components contains degenerated triangles");
+    }
+#endif
+  }
+  tree.insert(triangles.begin(), triangles.end());
+  tree.accelerate_distance_queries();
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::copyleft::cgal::submesh_aabb_tree<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>, CGAL::Epeck>(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&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epeck, CGAL::AABB_triangle_primitive<CGAL::Epeck, std::__1::vector<CGAL::Epeck::Triangle_3, std::__1::allocator<CGAL::Epeck::Triangle_3> >::iterator, CGAL::Boolean_tag<false> > > >&, std::__1::vector<CGAL::Epeck::Triangle_3, std::__1::allocator<CGAL::Epeck::Triangle_3> >&, std::__1::vector<bool, std::__1::allocator<bool> >&);
+#endif

include/igl/copyleft/cgal/submesh_aabb_tree.h

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+//
+#ifndef IGL_COPYLET_CGAL_SUBMESH_AABB_TREE_H
+#define IGL_COPYLET_CGAL_SUBMESH_AABB_TREE_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/intersections.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Build an AABB tree for a submesh indicated by a face selection list I
+      // of a full mesh (V,F)
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertices.
+      //   F  #F by 3 array of faces.
+      //   I  #I list of triangle indices to consider.
+      // Outputs:
+      //   tree  aabb containing triangles of (V,F(I,:))
+      //   triangles  #I list of cgal triangles
+      //   in_I  #F list of whether in submesh
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename Kernel>
+      IGL_INLINE void submesh_aabb_tree(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        CGAL::AABB_tree<
+          CGAL::AABB_traits<
+            Kernel, 
+            CGAL::AABB_triangle_primitive<
+              Kernel, typename std::vector<
+                typename Kernel::Triangle_3 >::iterator > > > & tree,
+        std::vector<typename Kernel::Triangle_3 > & triangles,
+        std::vector<bool> & in_I);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "submesh_aabb_tree.cpp"
+#endif
+
+#endif