add intersection of two AABB_tree and parallel build

AABB_tree/benchmark/AABB_tree/CMakeLists.txt

@@ -6,8 +6,10 @@ project(AABB_traits_benchmark)
 
 find_package(CGAL REQUIRED)
 
-create_single_source_cgal_program("test.cpp")
+create_single_source_cgal_program("old_bench_AABB_tree.cpp")
 create_single_source_cgal_program("tree_construction.cpp")
+create_single_source_cgal_program("knot_generation.cpp")
+create_single_source_cgal_program("tree_queries.cpp")
 
 # google benchmark
 find_package(benchmark QUIET)
@@ -17,3 +19,11 @@ if(benchmark_FOUND)
 else()
   message(STATUS "NOTICE: The benchmark 'tree_creation.cpp' requires the Google benchmark library, and will not be compiled.")
 endif()
+
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(tree_construction PRIVATE CGAL::TBB_support)
+else()
+  message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.")
+endif()

AABB_tree/benchmark/AABB_tree/knot_generation.cpp

@@ -0,0 +1,122 @@
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Vector = Kernel::Vector_3;
+using Triangle = std::array<int, 3>;
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = argc > 1 ? argv[1] : "knot.off";
+
+  const int p = argc > 2 ? std::atoi(argv[2]) : 1; // Nb of major turns
+  const int q = argc > 3 ? std::atoi(argv[3]) : 6; // Nb of knot turns
+
+  const double R = argc > 4 ? std::atof(argv[4]) : 2.0; // major radius
+  const double r = argc > 5 ? std::atof(argv[5]) : 0.9; // knot radius
+  const double tube_radius = argc > 6 ? std::atof(argv[6]) : 0.4;
+
+  const int tubular_segments = argc > 7 ? std::atoi(argv[7]) : 120;
+  const int radial_segments  = argc > 8 ? std::atoi(argv[8]) : 60;
+
+  const double exc = argc > 9 ? std::atof(argv[9]) : 1; // excentricity
+
+  std::vector<Point> vertices;
+  std::vector<Triangle> faces;
+
+  // Generate knot polyline
+  auto knot = [&](double t){
+    double cq = std::cos(q * t);
+    double sq = std::sin(q * t);
+    double factor = R + r * cq;
+    return Point(factor * std::cos(p * t),
+                 factor * std::sin(p * t),
+                 r * sq);
+  };
+  auto knot_derivative = [&](double t){
+    double cq = std::cos(q * t);
+    double sq = std::sin(q * t);
+    double cp = std::cos(p * t);
+    double sp = std::sin(p * t);
+
+    double factor = R + r*cq;
+    double dfactor = -r*q*sq;
+
+    return Vector(dfactor*cp - factor*p*sp,
+                  dfactor*sp + factor*p*cp,
+                  r*q*cq);
+  };
+
+  const double dt = 2.0 * M_PI / tubular_segments;
+  const double dr = 2.0 * M_PI / radial_segments;
+
+  // A normal to the knot axis
+  Vector N(1,0,0);
+  for(int i=0; i<tubular_segments; ++i)
+  {
+    double t = i * dt;
+
+    // knot point and its tangent
+    Point C  = knot(t);
+    Vector T = knot_derivative(t);
+    // Compute bases of the plane on C
+    Vector B =  CGAL::cross_product(N, T);
+    N = -CGAL::cross_product(B, T);
+    N /= CGAL::approximate_sqrt(N.squared_length());
+    B /= CGAL::approximate_sqrt(B.squared_length());
+
+    for(int j=0; j<radial_segments; ++j){
+      double phi = j * dr;
+      vertices.push_back(C + tube_radius*(exc*std::cos(phi)*N + 1/exc*std::sin(phi)*B));
+    }
+  }
+
+  auto idx = [&](int i, int j){
+    i %= tubular_segments;
+    j %= radial_segments;
+    return i * radial_segments + j;
+  };
+
+  for(int i=0; i<tubular_segments-1; ++i)
+    for(int j=0; j<radial_segments; ++j){
+      int v00 = idx(i  , j  );
+      int v10 = idx(i+1, j  );
+      int v01 = idx(i  , j+1);
+      int v11 = idx(i+1, j+1);
+
+      faces.push_back({v00,v10,v11});
+      faces.push_back({v00,v11,v01});
+    }
+
+  // To avoid a twist, we compute the offset that minimize the distance between a point at tubular index N-1 and tubular index 0
+  int n = tubular_segments-1;
+  int offset = 0;
+  double min = CGAL::squared_distance(vertices[idx(n, 0)], vertices[idx(0,0)]);
+  for(int j=1; j<radial_segments; ++j){
+    double sq = CGAL::squared_distance(vertices[idx(n, 0)], vertices[idx(0,j)]);
+    if(min >sq){
+      min = sq;
+      offset = j;
+    }
+  }
+  for(int j=0; j<radial_segments; ++j){
+    int v00 = idx(n, j  );
+    int v10 = idx(0, j+offset);
+    int v01 = idx(n, j+1);
+    int v11 = idx(0, j+offset+1);
+
+    faces.push_back({v00,v10,v11});
+    faces.push_back({v00,v11,v01});
+  }
+
+  CGAL::IO::write_polygon_soup(filename, vertices, faces);
+
+  return 0;
+}

AABB_tree/benchmark/AABB_tree/test.cpp → AABB_tree/benchmark/AABB_tree/old_bench_AABB_tree.cpp

@@ -8,6 +8,7 @@
 #include <CGAL/Polygon_mesh_processing/transform.h>
 #include <CGAL/Rigid_triangle_mesh_collision_detection.h>
 #include <CGAL/Side_of_triangle_mesh.h>
+#include <CGAL/Real_timer.h>
 
 #include <fstream>
 #include <sstream>
@@ -157,17 +158,16 @@ int main(int argc, const char** argv)
   std::string path = (argc>2)?argv[2]: CGAL::data_file_path("meshes/handle.off");
 
   std::cout<< k<<" steps in "<<path<<std::endl;
+  CGAL::Real_timer t;
   int nb_inter(0), nb_no_inter(0), nb_include(0),
       naive_inter(0), naive_no_inter(0), naive_include(0);
-  auto start = std::chrono::steady_clock::now();
-  naive_test(k, path, naive_inter, naive_no_inter, naive_include);
-  auto end = std::chrono::steady_clock::now();
+  t.start();
+      naive_test(k, path, naive_inter, naive_no_inter, naive_include);
   std::cout<<"Naive test           : "<<naive_inter<<" collisions, "<<naive_include<<" inclusions, "<<naive_no_inter<<" no collision, calculated in "
-          <<std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms." << std::endl;
-  start = std::chrono::steady_clock::now();
+          << t.time() << " ms." << std::endl;
+  t.stop(); t.reset(); t.start();
   test_no_collision(k, path,nb_inter, nb_no_inter, nb_include);
-  end = std::chrono::steady_clock::now();
   std::cout<<"With transform_traits: "<<nb_inter<<" collisions, "<<nb_include<<" inclusions, "<<nb_no_inter<<" no collision, calculated in "
-          <<std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms." << std::endl;
+          << t.time() << " ms." << std::endl;
     return 0;
 }

AABB_tree/benchmark/AABB_tree/tree_construction.cpp

@@ -1,12 +1,14 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/AABB_tree.h>
-#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_traits_3.h>
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 #include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
 
 #include <iostream>
 #include <fstream>
@@ -61,30 +63,52 @@ struct Compute_bbox {
   BBM bbm;
 };
 
-template <class K>
+template <typename ConcurrencyTag, class K>
 void run(std::string input)
 {
   typedef typename K::Point_3 Point_3;
   typedef CGAL::Surface_mesh<Point_3> Mesh;
   typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
-  typedef CGAL::AABB_traits<K, Primitive> Traits;
+  typedef CGAL::AABB_traits_3<K, Primitive> Traits;
   typedef CGAL::AABB_tree<Traits> Tree;
 
   Mesh tm;
-  std::ifstream(input) >> tm;
+  CGAL::IO::read_polygon_mesh(input, tm);
 
   {
   Tree tree(faces(tm).begin(), faces(tm).end(), tm);
-  CGAL::Timer time;
+  CGAL::Real_timer time;
   time.start();
-  tree.build();
-  time.stop();
+  tree.template build<ConcurrencyTag>();
   std::cout << "  build() time: " << time.time() << "\n";
+  tree.template accelerate_distance_queries<ConcurrencyTag>();
+  std::cout << "  build() + build kd-tree time: " << time.time() << "\n";
+  }
+
+  {
+  typedef CGAL::dynamic_face_property_t<CGAL::Bbox_3> Face_bbox_tag;
+  typedef typename boost::property_map<Mesh,  Face_bbox_tag >::type BboxMap;
+  typedef CGAL::AABB_traits_3<K, Primitive, BboxMap> BTraits;
+  typedef CGAL::AABB_tree<BTraits> BTree;
+
+  auto bb = get( Face_bbox_tag(), tm);
+  for(auto fd : faces(tm))
+      put(bb, fd, CGAL::Polygon_mesh_processing::face_bbox(fd, tm));
+
+  BTraits traits(bb);
+  BTree tree(traits);
+  tree.insert(faces(tm).begin(), faces(tm).end(), tm);
+  CGAL::Real_timer time;
+  time.start();
+  tree.template build<ConcurrencyTag>();
+  std::cout << "  build() with reference bbox time: " << time.time() << "\n";
+  tree.template accelerate_distance_queries<ConcurrencyTag>();
+  std::cout << "  build() with reference bbox + build kd-tree time: " << time.time() << "\n";
   }
 
   {
   Tree tree(faces(tm).begin(), faces(tm).end(), tm);
-  CGAL::Timer time;
+  CGAL::Real_timer time;
   time.start();
 
   typedef CGAL::Pointer_property_map<CGAL::Bbox_3>::type BBM;
@@ -108,17 +132,26 @@ void run(std::string input)
 
   Compute_bbox<BBM> compute_bbox(bbm);
   Split_primitives<RPM> split_primitives(rpm);
-  tree.custom_build(compute_bbox, split_primitives);
-  time.stop();
+  tree.template custom_build<ConcurrencyTag>(compute_bbox, split_primitives);
   std::cout << "  custom_build() time: " << time.time() << "\n";
+  tree.template accelerate_distance_queries<ConcurrencyTag>();
+  std::cout << "  custom_build() + build kd-tree time: " << time.time() << "\n";
   }
 }
 
 int main(int, char** argv)
 {
+  std::cout << "Build with Cartesian\n";
+  run<CGAL::Parallel_tag, CGAL::Simple_cartesian<double>>(argv[1]);
   std::cout << "Build with Epick\n";
-  run<CGAL::Epick>(argv[1]);
+  run<CGAL::Parallel_tag, CGAL::Epick>(argv[1]);
   std::cout << "Build with Epeck\n";
-  run<CGAL::Epeck>(argv[1]);
+  run<CGAL::Parallel_tag, CGAL::Epeck>(argv[1]);
+  std::cout << "Build with Cartesian (Sequential) \n";
+  run<CGAL::Sequential_tag, CGAL::Simple_cartesian<double>>(argv[1]);
+  std::cout << "Build with Epick (Sequential) \n";
+  run<CGAL::Sequential_tag, CGAL::Epick>(argv[1]);
+  std::cout << "Build with Epeck (Sequential) \n";
+  run<CGAL::Sequential_tag, CGAL::Epeck>(argv[1]);
   return EXIT_SUCCESS;
 }

AABB_tree/benchmark/AABB_tree/tree_creation.cpp

@@ -55,7 +55,7 @@ BENCHMARK(BM_Intersections);
 int main(int argc, char** argv)
 {
   std::string default_file = CGAL::data_file_path("meshes/handle.off");
-  std::strint filename = argc > 2? argv[2] : default_file;
+  std::string filename = argc > 2? argv[2] : default_file;
 
   {
     std::ifstream input(filename);