diff --git a/AABB_tree/benchmark/AABB_tree/CMakeLists.txt b/AABB_tree/benchmark/AABB_tree/CMakeLists.txt
index bfcbb5454fa1..4e2b57075b1a 100644
--- a/AABB_tree/benchmark/AABB_tree/CMakeLists.txt
+++ b/AABB_tree/benchmark/AABB_tree/CMakeLists.txt
@@ -6,8 +6,10 @@ project(AABB_traits_benchmark)
 
 find_package(CGAL REQUIRED OPTIONAL_COMPONENTS Core)
 
-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()
diff --git a/AABB_tree/benchmark/AABB_tree/knot_generation.cpp b/AABB_tree/benchmark/AABB_tree/knot_generation.cpp
new file mode 100644
index 000000000000..aa19db304e44
--- /dev/null
+++ b/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;
+}
\ No newline at end of file
diff --git a/AABB_tree/benchmark/AABB_tree/test.cpp b/AABB_tree/benchmark/AABB_tree/old_bench_AABB_tree.cpp
similarity index 92%
rename from AABB_tree/benchmark/AABB_tree/test.cpp
rename to AABB_tree/benchmark/AABB_tree/old_bench_AABB_tree.cpp
index 5ef17a2cf317..c0a6c757975f 100644
--- a/AABB_tree/benchmark/AABB_tree/test.cpp
+++ b/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;
 }
diff --git a/AABB_tree/benchmark/AABB_tree/tree_construction.cpp b/AABB_tree/benchmark/AABB_tree/tree_construction.cpp
index 3ce2948b0506..2ab8d413a4f3 100644
--- a/AABB_tree/benchmark/AABB_tree/tree_construction.cpp
+++ b/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;
 }
diff --git a/AABB_tree/benchmark/AABB_tree/tree_creation.cpp b/AABB_tree/benchmark/AABB_tree/tree_creation.cpp
index c5bd5396d513..fe074cb5e3cf 100644
--- a/AABB_tree/benchmark/AABB_tree/tree_creation.cpp
+++ b/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);
diff --git a/AABB_tree/benchmark/AABB_tree/tree_queries.cpp b/AABB_tree/benchmark/AABB_tree/tree_queries.cpp
new file mode 100644
index 000000000000..b27e44b714f3
--- /dev/null
+++ b/AABB_tree/benchmark/AABB_tree/tree_queries.cpp
@@ -0,0 +1,228 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Random.h>
+#include <CGAL/Real_timer.h>
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+using K = CGAL::Simple_cartesian<double>;
+
+using Point   = K::Point_3;
+using Vector  = K::Vector_3;
+using Segment = K::Segment_3;
+using Ray     = K::Ray_3;
+using Line    = K::Line_3;
+using Plane   = K::Plane_3;
+
+using Mesh = CGAL::Surface_mesh<Point>;
+
+using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
+using Traits    = CGAL::AABB_traits_3<K, Primitive>;
+using Tree      = CGAL::AABB_tree<Traits>;
+
+CGAL::Random rng;
+
+template<class Kernel = K>
+typename Kernel::Point_3 random_point(const CGAL::Bbox_3& bb){
+  return typename Kernel::Point_3(rng.get_double(bb.xmin(), bb.xmax()),
+                                  rng.get_double(bb.ymin(), bb.ymax()),
+                                  rng.get_double(bb.zmin(), bb.zmax()));
+}
+
+Vector random_vector(){
+  Vector v;
+  do{
+    v = Vector(rng.get_double(-1., 1.),
+               rng.get_double(-1., 1.),
+               rng.get_double(-1., 1.));
+  } while(v.squared_length() < 1e-12);
+  return v;
+}
+
+template<class QueryGenerator, class Function>
+double benchmark(QueryGenerator gen,
+                 Function f,
+                 std::size_t nb_queries = 500000)
+{
+  CGAL::Real_timer timer;
+  timer.start();
+  for(std::size_t i=0; i<nb_queries; ++i)
+    f(gen());
+  timer.stop();
+  return double(nb_queries) / timer.time();
+}
+
+void benchmark_queries(const Mesh& mesh)
+{
+  Tree tree(faces(mesh).first, faces(mesh).second, mesh);
+  tree.build();
+
+  CGAL::Bbox_3 bb = PMP::bbox(mesh);
+  constexpr std::size_t N = 100000;
+  // Warm-up
+  for(int i=0;i<10000;++i)
+    tree.do_intersect(Segment(random_point(bb), random_point(bb)));
+  std::cout << std::fixed << std::setprecision(0);
+
+  auto print = [](const std::string& name, double value)  {
+    std::cout << std::setw(40) << std::left
+              << name
+              << value
+              << " queries/s\n";
+  };
+
+#define BENCHMARK(QueryType, Generator)                                      \
+  {                                                                          \
+    std::cout << "\n=== " #QueryType " ===\n";                               \
+                                                                             \
+    print("do_intersect",                                                    \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      { tree.do_intersect(q); }, N));                                        \
+                                                                             \
+    print("any_intersected_primitive",                                       \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      { tree.any_intersected_primitive(q); }, N));                           \
+                                                                             \
+    print("any_intersection",                                                \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      { tree.any_intersection(q); }, N));                                    \
+                                                                             \
+    print("number_of_intersected_primitives",                                \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      { tree.number_of_intersected_primitives(q); }, N));                    \
+                                                                             \
+    print("all_intersected_primitives",                                      \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      {                                                                      \
+        std::vector<Tree::Primitive_id> ids;                                 \
+        tree.all_intersected_primitives(q, std::back_inserter(ids));          \
+      }, N));                                                                \
+                                                                             \
+    print("all_intersections",                                               \
+      benchmark(Generator,                                                   \
+      [&](const QueryType& q)                                                \
+      {                                                                      \
+        using Intersection =                                                 \
+          typename Tree::template Intersection_and_primitive_id<QueryType>::Type; \
+                                                                             \
+        std::vector<Intersection> out;                                       \
+        tree.all_intersections(q, std::back_inserter(out));                  \
+      }, N));                                                                \
+  }
+
+  BENCHMARK(Segment, [&](){ return Segment(random_point(bb), random_point(bb)); });
+  BENCHMARK(Ray, [&](){ return Ray(random_point(bb), random_vector()); });
+  BENCHMARK(Line, [&](){ return Line(random_point(bb), random_point(bb)); });
+  BENCHMARK(Plane, [&](){ return Plane(random_point(bb), random_vector()); });
+#undef BENCHMARK
+}
+
+template<class Kernel>
+void benchmark_kernel(const std::string &filename, const std::string &kernel_name)
+{
+  using Segment = typename Kernel::Segment_3;
+
+  using Mesh = CGAL::Surface_mesh<typename Kernel::Point_3>;
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
+  using Traits    = CGAL::AABB_traits_3<Kernel, Primitive>;
+  using Tree      = CGAL::AABB_tree<Traits>;
+  const int N = 100000;
+
+  Mesh mesh;
+  CGAL::IO::read_polygon_mesh(filename, mesh);
+  Tree tree(faces(mesh).first, faces(mesh).second, mesh);
+  tree.build();
+  CGAL::Bbox_3 bb = PMP::bbox(mesh);
+
+  CGAL::Real_timer t;
+  t.start();
+  for(std::size_t i=0;i<N;++i){
+      Segment s(random_point<Kernel>(bb), random_point<Kernel>(bb));
+      using Intersection = typename Tree::template Intersection_and_primitive_id<Segment>::Type;
+      std::vector<Intersection> out;
+      tree.all_intersections(s, std::back_inserter(out));
+  }
+  t.stop();
+
+  std::cout << std::setw(40)
+            << kernel_name
+            << N/t.time()
+            << '\n';
+}
+
+void benchmark_distances(const Mesh& mesh, std::size_t N = 300000)
+{
+  Tree tree(faces(mesh).first, faces(mesh).second, mesh);
+  tree.build();
+  tree.accelerate_distance_queries();
+
+  const CGAL::Bbox_3 bb = PMP::bbox(mesh);
+
+  std::vector<Point> queries;
+  queries.reserve(N);
+  for(std::size_t i=0; i<N; ++i)
+    queries.push_back(random_point(bb));
+
+  auto benchmark = [&](const std::string& name, auto&& f)  {
+    CGAL::Real_timer timer;
+    timer.start();
+    for(const Point& p : queries)
+      f(p);
+    timer.stop();
+
+    std::cout << std::setw(32) << std::left
+              << name
+              << std::fixed << std::setprecision(0)
+              << N / timer.time()
+              << " queries/s\n";
+  };
+
+  benchmark("closest_point()", [&](const Point& p){ tree.closest_point(p); });
+  benchmark("squared_distance()", [&](const Point& p){ tree.squared_distance(p); });
+  benchmark("closest_point_and_primitive()", [&](const Point& p){ tree.closest_point_and_primitive(p); });
+}
+
+int main(int argc, char** argv)
+{
+  if(argc != 2){
+    std::cerr << "Usage: " << argv[0] << " mesh.off\n";
+    return EXIT_FAILURE;
+  }
+
+  Mesh mesh;
+  std::string filename = argv[1];
+
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh))  {
+    std::cerr << "Cannot read mesh\n";
+    return EXIT_FAILURE;
+  }
+  // benchmark_queries(mesh);
+  // benchmark_kernel<CGAL::Simple_cartesian<double>>(filename, "Simple_cartesian<double>");
+  // benchmark_kernel<CGAL::Simple_cartesian<float>>(filename, "Simple_cartesian<float>");
+  // benchmark_kernel<CGAL::Cartesian<double>>(filename, "Cartesian<double>");
+  // benchmark_kernel<CGAL::Cartesian<float>>(filename, "Cartesian<float>>");
+  // benchmark_kernel<CGAL::Epick>(filename, "Epick");
+  benchmark_distances(mesh);
+
+  return EXIT_SUCCESS;
+}
\ No newline at end of file
diff --git a/AABB_tree/doc/AABB_tree/Doxyfile.in b/AABB_tree/doc/AABB_tree/Doxyfile.in
index fb005db4f918..5ea1773a1715 100644
--- a/AABB_tree/doc/AABB_tree/Doxyfile.in
+++ b/AABB_tree/doc/AABB_tree/Doxyfile.in
@@ -5,3 +5,5 @@ PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Fast Intersection and Distance Comput
 EXTRACT_ALL                =  false
 HIDE_UNDOC_MEMBERS         =  true
 HIDE_UNDOC_CLASSES         =  true
+
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/AABB_trees/intersection.h
diff --git a/AABB_tree/doc/AABB_tree/PackageDescription.txt b/AABB_tree/doc/AABB_tree/PackageDescription.txt
index e16cfc867d40..3f6361165814 100644
--- a/AABB_tree/doc/AABB_tree/PackageDescription.txt
+++ b/AABB_tree/doc/AABB_tree/PackageDescription.txt
@@ -50,4 +50,8 @@
 - `CGAL::AABB_primitive<Id,ObjectPropertyMap,PointPropertyMap,ExternalPropertyMaps,CacheDatum>`
 - `CGAL::AABB_halfedge_graph_segment_primitive<HalfedgeGraph,Vpm,OneHalfedgeGraphPerTree,CacheDatum>`
 - `CGAL::AABB_face_graph_triangle_primitive<FaceGraph,Vpm,OneFaceGraphPerTree,CacheDatum>`
+
+\cgalCRPSection{Functions}
+- `CGAL::AABB_trees::do_intersect()`
+- `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()`
 */
diff --git a/AABB_tree/doc/AABB_tree/aabb_tree.txt b/AABB_tree/doc/AABB_tree/aabb_tree.txt
index fbdc0b890c32..f4eba8f292ea 100644
--- a/AABB_tree/doc/AABB_tree/aabb_tree.txt
+++ b/AABB_tree/doc/AABB_tree/aabb_tree.txt
@@ -23,10 +23,10 @@ triangles, or plane objects (planes, triangles) against sets of
 segments. An example of a distance query consists of finding the
 closest point from a point query to a set of triangles.
 
-Note that this component is not suited to the problem of finding all
-intersecting pairs of objects. We refer to the component
+This component solves the problem of finding all
+intersecting pairs of objects and achieves better perfomance than
 \ref chapterBoxIntersection "Intersecting Sequences of dD Iso-oriented Boxes"
-which can find all intersecting pairs of iso-oriented boxes.
+when using parallel execution.
 
 The AABB tree data structure takes as input an iterator range of
 geometric data, which is then converted into primitives. From these
@@ -109,6 +109,18 @@ For example if one is using `CGAL::AABB_traits` with a Kernel from \cgal,
 having degenerate triangles or segments in the AABB-tree will result in an undefined
 behavior or a crash.
 
+\section two_aabb_tree_interface Operations between Two AABB Trees
+
+Given two AABB trees, the function `CGAL::AABB_trees::do_intersect()` determines whether any
+primitive of the first tree intersects a primitive of the second tree.
+To enumerate all intersecting pairs of primitives, use `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()`.
+
+These functions provide an alternative to `CGAL::box_intersection_d()`. In practice, `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()` generally provides better performance,
+particularly when using a parallel execution tag (see \ref aabb_tree_perf).
+
+Moreover, if several calls to `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()` are performed using the same trees, the trees need to be constructed only once.
+Consequently, the performance advantage over `CGAL::box_intersection_d()` becomes even greater.
+
 \section aabb_tree_examples Examples
 
 \subsection aabb_tree_examples_1 Tree of Triangles, for Intersection and Distance Queries
@@ -200,22 +212,29 @@ to custom points. In \ref AABB_tree/AABB_custom_triangle_soup_example.cpp "AABB_
 triangles are stored into a single array as to form a triangle
 soup. The primitive internally uses a `boost::iterator_adaptor`
 as to provide the three functions `AABBPrimitive::id()`, `AABBPrimitive::datum()`,
-and `AABBPrimitive::reference_point()`) required by the primitive concept. In
+and `AABBPrimitive::reference_point()` required by the primitive concept. In
 \ref AABB_tree/AABB_custom_indexed_triangle_set_example.cpp "AABB_custom_indexed_triangle_set_example.cpp" the input is an indexed
 triangle set stored through two arrays: one array of points and one
 array of indices which refer to the point array. Here also the
 primitive internally uses a `boost::iterator_adaptor`.
 
+\subsection aabb_tree_examples_9 Intersection of Two Trees
+
+In the following example, we compute all the intersections between two tetrahedra. The tetrahedra
+are stored as a vector of `Triangle_3`. We first compute if the tetrahedra do intersect and then we
+compute the indices of the intersecting triangles.
+
+\cgalExample{AABB_tree/AABB_intersection_example.cpp}
+
 \section aabb_tree_perf Performances
 
 We provide some performance numbers for the case where the AABB tree
 contains a set of polyhedron triangle facets. We measure the tree
 construction time, the memory occupancy and the number of queries per
 second for a variety of intersection and distance queries. The machine
-used is a PC running Windows XP64 with an Intel CPU Core2 Extreme
-clocked at 3.06 GHz with 4GB of RAM. By default, the kernel used is
-`Simple_cartesian<double>` (the fastest in our experiments). The
-program has been compiled with Visual C++ 2005 compiler with the O2
+used is a PC running Linux Mint with an Intel CPU Core Ultra 7 155H with 16 cores and 32GB of RAM.
+By default, the kernel used is `Simple_cartesian<double>` (the fastest in our experiments). The
+program has been compiled with gcc 13.3 2023 compiler with the O3
 option which maximizes speed.
 
 \subsection aabb_tree_perf_cons Construction
@@ -223,22 +242,20 @@ option which maximizes speed.
 The surface triangle mesh chosen for benchmarking the tree
 construction is the knot model (14,400 triangles) depicted by
 \cgalFigureRef{figAABB-tree-bench}. We measure the tree construction time (both
-AABB tree alone and AABB tree with internal KD-tree) for this model as
-well as for three denser versions subdivided through the Loop
-subdivision scheme which multiplies the number of triangles by four.
+AABB tree alone and AABB tree with internal KD-tree) for four denser versions if the knot model.
 
 
-|  Number of triangles | Construction (in ms) | Construction with internal KD-tree (in ms) |
-| ----: | ----:     | ----:     |
-   14,400 |     156 |       157 |
-   57,600 |     328 |       328 |
-  230,400 |   1,141 |     1,437 |
-  921,600 |   4,813 |     5,953 |
+|  Number of triangles | Construction (in ms) (Parallel) | (Sequential) | Construction with internal KD-tree (in ms) (Parallel) | (Sequential) |
+| ----: | ----: | ----: | ----: | ----: |
+  230,400 | 20 | 36 | 48 | 66 |
+  921,600 | 72 | 187 | 155 | 306 |
+  3,686,400 | 292 | 816 | 729 | 1426 |
+  14,745,600 | 1,425 | 4,196 | 3,028 | 7,016 |
 
 \subsection aabb_tree_perf_mem Memory
 
 When using the polyhedron triangle facet primitive (defined in
-`AABB_face_graph_triangle_primitive.h`) the AABB tree occupies
+`CGAL/AABB_face_graph_triangle_primitive.h`) the AABB tree occupies
 approximately 61 bytes per primitive (without constructing the
 internal KD-tree). It increases to approximately 150 bytes per
 primitive when constructing the internal KD-tree with one reference
@@ -269,7 +286,7 @@ takes an iterator range as input.
 \subsection aabb_tree_perf_inter Intersections
 
 The following table measures the number of intersection queries per
-second on the 14,400 triangle version of the knot mesh model for ray,
+second on one core on the 14,400 triangle version of the knot mesh model for ray,
 line, segment and plane queries. Each ray query is generated by
 choosing a random source point within the mesh bounding box and a
 random vector. A line or segment query is generated by choosing two
@@ -282,45 +299,38 @@ the intersection functions which enumerate all intersections.
 
 |  Function                              | Segment |     Ray |    Line |   Plane |
 | :---- | ----: | ----: | -: | -: |
-|   AABB_tree::do_intersect()                       | 187,868 | 185,649 | 206,096 | 377,969 |
-|   AABB_tree::any_intersected_primitive()         | 190,684 | 190,027 | 208,941 | 360,337 |
-|   AABB_tree::any_intersection()                   | 147,468 | 143,230 | 148,235 | 229,336 |
-|   AABB_tree::number_of_intersected_primitives() |  64,389 |  52,943 |  54,559 |   7,906 |
-|   AABB_tree::all_intersected_primitives()        |  65,553 |  54,838 |  53,183 |   5,693 |
-|   AABB_tree::all_intersections()                  |  46,507 |  38,471 |  36,374 |   2,644 |
-
+|   AABB_tree::do_intersect()                       | 513,505 | 788,246 | 773,947 | 1,863,721  |
+|   AABB_tree::any_intersected_primitive()         | 521,651 | 852,929 | 831,919 | 2,125,582 |
+|   AABB_tree::any_intersection()                   | 504,836 | 835,164 | 821,369 | 1,730,017 |
+|   AABB_tree::number_of_intersected_primitives() |  186,102 |  397,069 |  210,052 |   49,941 |
+|   AABB_tree::all_intersected_primitives()        |  178,769 |  397,458 |  192,482 |   49,324 |
+|   AABB_tree::all_intersections()                  |  178,536 |  352,541 |  193,700 |   21,745 |
 
 Curve of \cgalFigureRef{figAABB-tree-bench} plots the number of queries per second
 (here the `AABB_tree::all_intersections()` function with random segment
 queries) against the number of input triangles for the knot triangle
 surface mesh.
 
-
-\cgalFigureBegin{figAABB-tree-bench,bench.png}
-Number of queries per second against number of triangles for the knot model with 14K (shown), 57K, 230K and 921K triangles. We call the `all_intersections()` function with segment queries randomly chosen within the bounding box.
+\cgalFigureBegin{figAABB-tree-bench,bench.png, knot.png}
+Number of queries per second against number of triangles for the knot model with 14K (shown) to 15M triangles. We call the `all_intersections()` function with segment queries randomly chosen within the bounding box.
 \cgalFigureEnd
 
 The following table measures the number of `AABB_tree::all_intersections()`
 queries per second against several kernels. We use the 14,400 triangle
-version of the knot mesh model for random segment queries. Note how
-the `Simple_cartesian` kernel is substantially faster than the
-`Cartesian` kernel.
+version of the knot mesh model for random segment queries.
 
 |  Kernel | Queries/s (all_intersections() with segment queries)|
 | :---- | ----: |
-|`Simple_cartesian<double>`               | 46,507 |
-|`Simple_cartesian<float>`                | 43,187 |
-|`Cartesian<double>`                      |  5,335 |
-|`Cartesian<float>`                       |  5,522 |
-|`Exact_predicates_inexact_constructions_kernel` | 18,411 |
-
-
+|`Simple_cartesian<double>`               | 201,618 |
+|`Simple_cartesian<float>`                | 193,414 |
+|`Cartesian<double>`                      | 164,784 |
+|`Cartesian<float>`                       | 159,686 |
+|`Exact_predicates_inexact_constructions_kernel` | 119,590 |
 
 \subsection aabb_tree_perf_dist Distances
 
 The surface triangle mesh chosen for benchmarking distances is again
-the knot model in four increasing resolutions obtained through Loop
-subdivision. In the following table we first measure the tree
+the knot model in six increasing resolutions. In the following table we first measure the tree
 construction time (which includes the construction of the internal
 KD-tree data structure used to accelerate the distance queries by up
 to one order of magnitude in our experiments). We then measure the
@@ -329,12 +339,36 @@ number of queries per second for the three types distance queries
 `AABB_tree::closest_point_and_primitive()`) from point queries randomly chosen
 inside the bounding box.
 
-|  Nb triangles | Construction (ms) | Closest_point() | Squared_distance() | Closest_point_and_primitive() |
-| ----:    | ----:     | ----: | ----: | -: |
-|   14,400 |   157.000 | 45,132 | 45,626 | 45,770 |
-|   57,600 |   328.000 | 21,589 | 21,312 | 21,137 |
-|  230,400 |     1.437 | 11,063 | 10,962 | 11,086 |
-|  921,600 |     5.953 |  5,636 |  5,722 |  5,703 |
+|  Nb triangles | Closest_point() | Squared_distance() | Closest_point_and_primitive() |
+| ----:    | ----: | ----: | -: |
+|   14,400 | 196,043 | 196,356 | 196,708 |
+|   57,600 | 123,388 | 123,844 | 123,567 |
+|  230,400 | 54,480 | 53,235 | 54,329 |
+|  921,600 | 25,761 | 26,890 | 25,386 |
+|  3,686,400 | 13,696 | 13,389 | 13,440 |
+|  14,745,600 | 7,158 | 6,913| 7,010 |
+
+\subsection aabb_two_trees_perf Intersection of Two AABB Trees
+
+We test the intersection of two AABB trees using 3 models: iphigenia and the Nefertiti that are two common models of mesh processing and twisted knot to maximize the number of intersection. For each model,
+we test their intersection with a translated version of themselves (Note: We do not exploit the fact that it is the same model, the tree is built two times).
+We compare here `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()` and `CGAL::box_intersection_d()`.
+
+We exhibit the performance of the intersection of two AABB trees using three models: Iphigenia and Nefertiti \cgalFigureRef{figAABB-tree-bench}, which are widely used benchmark meshes in geometry processing, and a twisted knot \cgalFigureRef{figAABB-tree-bench}, chosen to maximize the number of intersections.
+
+For each model, we compute the intersections between the model and a translated copy of itself. Although both trees are built from the same mesh, we intentionally construct the AABB tree twice and do not exploit this fact.
+We compare the performance of `CGAL::AABB_trees::all_pairs_of_intersecting_primitives()` with that of `CGAL::box_intersection_d()`.
+
+\cgalFigureBegin{figAABB-tree-model,iphigenia.png, nefertiti.png, twisted_knot.png}
+The three input models used to illustrate AABB trees intersection performance. Iphigenia (left), Nefertiti (center), Twisted knot (right).
+\cgalFigureEnd
+
+| Models | Nb Triangles | Nb Intersections | **AABB** | Tree (ms) | | | **Box** d || **AABB** | Tree Seq. | (ms) | | **Box** d Seq. |
+| ----: | ----: | ----: | ----: | ----: | ----: | ----: | ----: || ----: | ----: | ----: | ----: | ----: |
+|  |  |  | Build |  Do intersect | All pairs | All pairs + Build |  || Build | Do intersect | All pairs | All pairs + Build |  |
+| Iphigenia × Iphigenia Tr. | 703,512 × 703,512 | 47,677 | 91 | <0.1 | 37 | 128 | 424 || 318 | <0.1 | 169 | 488 | 868 |
+| Nefertiti × Nefertiti Tr. | 2,018,232 × 2,018,232 | 14,246 | 261 | <0.1 | 10 | 271 | 1,006 ||1,087 | <0.1 | 46 | 1,134 | 2,348 |
+| Twisted knot × Twisted knot Rot. | 600,000 × 600,000 | 240,144 | 154 | <0.1 | 57 | 211 | 850 || 541 | <0.1 | 440 | 981 | 2,993 |
 
 \subsection aabb_tree_perf_summary Summary
 
@@ -450,6 +484,8 @@ thread-safe queries, introduction of shared data stored in the traits
 for lighter primitive types, ...
 In 2024, the package was made compatible with 2D and 3D primitives by
 Andreas Fabri, Sébastien Loriot, and Sven Oesau.
+In 2026, the package was extended to support the intersection of two AABB trees.
+This work also introduced performance optimizations and parallelization, contributed by Leo Valque.
 
 */
 
diff --git a/AABB_tree/doc/AABB_tree/examples.txt b/AABB_tree/doc/AABB_tree/examples.txt
index d986357796cd..7add888c9321 100644
--- a/AABB_tree/doc/AABB_tree/examples.txt
+++ b/AABB_tree/doc/AABB_tree/examples.txt
@@ -16,4 +16,5 @@
 \example AABB_tree/AABB_triangle_3_example.cpp
 \example AABB_tree/AABB_halfedge_graph_edge_example.cpp
 \example AABB_tree/AABB_face_graph_triangle_example.cpp
+\example AABB_tree/AABB_intersection_example.cpp
 */
diff --git a/AABB_tree/doc/AABB_tree/fig/bench.png b/AABB_tree/doc/AABB_tree/fig/bench.png
index b64ca29deef7..5d2fa041d20e 100644
Binary files a/AABB_tree/doc/AABB_tree/fig/bench.png and b/AABB_tree/doc/AABB_tree/fig/bench.png differ
diff --git a/AABB_tree/doc/AABB_tree/fig/iphigenia.png b/AABB_tree/doc/AABB_tree/fig/iphigenia.png
new file mode 100644
index 000000000000..589524a1ca84
Binary files /dev/null and b/AABB_tree/doc/AABB_tree/fig/iphigenia.png differ
diff --git a/AABB_tree/doc/AABB_tree/fig/knot.png b/AABB_tree/doc/AABB_tree/fig/knot.png
new file mode 100644
index 000000000000..215a63ff7f4b
Binary files /dev/null and b/AABB_tree/doc/AABB_tree/fig/knot.png differ
diff --git a/AABB_tree/doc/AABB_tree/fig/nefertiti.png b/AABB_tree/doc/AABB_tree/fig/nefertiti.png
new file mode 100644
index 000000000000..430d7362df8f
Binary files /dev/null and b/AABB_tree/doc/AABB_tree/fig/nefertiti.png differ
diff --git a/AABB_tree/doc/AABB_tree/fig/twisted_knot.png b/AABB_tree/doc/AABB_tree/fig/twisted_knot.png
new file mode 100644
index 000000000000..541dd108eb4b
Binary files /dev/null and b/AABB_tree/doc/AABB_tree/fig/twisted_knot.png differ
diff --git a/AABB_tree/examples/AABB_tree/AABB_intersection_example.cpp b/AABB_tree/examples/AABB_tree/AABB_intersection_example.cpp
new file mode 100644
index 000000000000..b31b3a2fdcd8
--- /dev/null
+++ b/AABB_tree/examples/AABB_tree/AABB_intersection_example.cpp
@@ -0,0 +1,59 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/AABB_triangle_primitive_3.h>
+#include <CGAL/AABB_trees/intersection.h>
+
+#include <iostream>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::Point_3 Point;
+typedef K::Triangle_3 Triangle;
+
+typedef std::vector<Triangle>::iterator Iterator;
+typedef CGAL::AABB_triangle_primitive_3<K, Iterator> Primitive;
+typedef CGAL::AABB_traits_3<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+int main()
+{
+    Point a(1.0, 0.0, 0.0);
+    Point b(0.0, 1.0, 0.0);
+    Point c(0.0, 0.0, 1.0);
+    Point d(0.0, 0.0, 0.0);
+
+    Point e(1.2, 0.2, 0.2);
+    Point f(0.2, 1.2, 0.2);
+    Point g(0.2, 0.2, 1.2);
+    Point h(0.2, 0.2, 0.2);
+
+    std::vector<Triangle> tetra1;
+    tetra1.push_back(Triangle(a,b,c));
+    tetra1.push_back(Triangle(a,b,d));
+    tetra1.push_back(Triangle(a,d,c));
+    tetra1.push_back(Triangle(b,c,d));
+
+    std::vector<Triangle> tetra2;
+    tetra2.push_back(Triangle(e,f,g));
+    tetra2.push_back(Triangle(e,f,h));
+    tetra2.push_back(Triangle(e,h,g));
+    tetra2.push_back(Triangle(f,g,h));
+
+    // constructs AABB tree
+    Tree tree1(tetra1.begin(),tetra1.end());
+    Tree tree2(tetra2.begin(),tetra2.end());
+
+    // do intersect
+    std::cout << "Tetrahedra do intersect: " << CGAL::AABB_trees::do_intersect(tree1, tree2) << std::endl;
+
+    std::vector< std::pair<Iterator, Iterator> > intersections;
+    CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(intersections));
+
+    for(auto [id1, id2]: intersections)
+      std::cout << "ids: " << std::distance(tetra1.begin(), id1) << " " << std::distance(tetra2.begin(), id2) << std::endl;
+
+    return EXIT_SUCCESS;
+}
diff --git a/AABB_tree/include/CGAL/AABB_indexed_triangle_primitive_3.h b/AABB_tree/include/CGAL/AABB_indexed_triangle_primitive_3.h
new file mode 100644
index 000000000000..aa5d55423e60
--- /dev/null
+++ b/AABB_tree/include/CGAL/AABB_indexed_triangle_primitive_3.h
@@ -0,0 +1,134 @@
+// Copyright (c) 2026 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Leo Valque
+//
+
+
+#ifndef CGAL_AABB_INDEXED_TRIANGLE_PRIMITIVE_3_H_
+#define CGAL_AABB_INDEXED_TRIANGLE_PRIMITIVE_3_H_
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/AABB_primitive.h>
+#include <iterator>
+
+namespace CGAL {
+
+namespace internal {
+
+template <class GeomTraits, class Iterator, class PointIterator, class PointMap>
+struct Triangle_3_from_index_range_iterator_property_map {
+  //classical typedefs
+  typedef Iterator key_type;
+  typedef typename GeomTraits::Triangle_3 value_type;
+  typedef typename GeomTraits::Triangle_3 reference;
+
+  typedef boost::readable_property_map_tag category;
+  typedef Triangle_3_from_index_range_iterator_property_map<GeomTraits, Iterator, PointIterator, PointMap> Self;
+
+  Triangle_3_from_index_range_iterator_property_map() {}
+  Triangle_3_from_index_range_iterator_property_map(PointIterator b, PointMap& pmap) : begin(b), pmap(pmap) {}
+
+  inline friend value_type
+  get(Self s, key_type it)
+  {
+    return typename GeomTraits::Construct_triangle_3()(get(s.pmap, s.begin[(*it)[0]]), get(s.pmap, s.begin[(*it)[1]]), get(s.pmap, s.begin[(*it)[2]]));
+  }
+
+  PointIterator begin;
+  PointMap pmap;
+};
+
+template <class GeomTraits, class Iterator, class PointIterator, class PointMap>
+struct Point_from_indexed_triangle_3_iterator_property_map {
+  //classical typedefs
+  typedef Iterator key_type;
+  typedef typename PointMap::value_type value_type;
+  typedef const value_type reference;
+
+  typedef boost::readable_property_map_tag category;
+  typedef Point_from_indexed_triangle_3_iterator_property_map<GeomTraits, Iterator, PointIterator, PointMap> Self;
+
+  Point_from_indexed_triangle_3_iterator_property_map() {}
+  Point_from_indexed_triangle_3_iterator_property_map(PointIterator b, PointMap &pmap) : begin(b), pmap(pmap) {}
+
+  inline friend reference
+  get(Self s, key_type it)
+  {
+    return get(s.pmap, s.begin[((*it)[0])]);
+  }
+
+  PointIterator begin;
+  PointMap pmap;
+};
+}//namespace internal
+
+
+/*!
+ * \ingroup PkgAABBTreeRef
+ * Primitive type that uses as identifier an iterator with a range of three indices as `value_type`.
+ * The iterator from which the primitive is built should not be invalided
+ * while the AABB tree holding the primitive is in use.
+ *
+ * \cgalModels{AABBPrimitive}
+ *
+ * \tparam GeomTraits is a traits class providing the nested type `Point_3` and `Triangle_3`.
+ *         It also provides the functor `Construct_triangle_3` that has an operator taking three `Point_3` as
+ *         parameters and returns a `Triangle_3`
+ * \tparam IndexIterator is a model of `ForwardIterator` with its value type being a `RandomAccessRange` of size 3 with an index type as `value_type`, e.g., `uint8_t`, `uint16_t` or int.
+ * \tparam PointRange is a model of `RandomAccessRange`. Its value type needs to be compatible to `PointMap` or `Point_3` in the default case.
+ * \tparam CacheDatum is either `CGAL::Tag_true` or `CGAL::Tag_false`. In the former case,
+ *         the datum is stored in the primitive, while in the latter it is
+ *         constructed on the fly to reduce the memory footprint.
+ *         The default is `CGAL::Tag_false` (datum is not stored).
+ * \tparam PointMap is a model of `ReadablePropertyMap` with its key type being the value type of `PointRange` and the value type being a `Point_3`.
+ *         The default is \link Identity_property_map `CGAL::Identity_property_map`\endlink<PointRange::value_type>.
+ *
+ * \sa `AABBPrimitive`
+ * \sa `AABB_primitive<Id,ObjectPropertyMap,PointPropertyMapPolyhedron,ExternalPropertyMaps,CacheDatum>`
+ * \sa `AABB_segment_primitive_2<GeomTraits,Iterator,CacheDatum>`
+ * \sa `AABB_triangle_primitive_2<GeomTraits,Iterator,CacheDatum>`
+ * \sa `AABB_triangle_primitive_3<GeomTraits,Iterator,CacheDatum>`
+ */
+template < class GeomTraits,
+           class IndexIterator,
+           class PointRange,
+           class CacheDatum = Tag_false,
+           class PointMap = Identity_property_map<typename PointRange::value_type>>
+class AABB_indexed_triangle_primitive_3
+#ifndef DOXYGEN_RUNNING
+  : public AABB_primitive<  IndexIterator,
+                            internal::Triangle_3_from_index_range_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>,
+                            internal::Point_from_indexed_triangle_3_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>,
+                            Tag_true,
+                            CacheDatum >
+#endif
+{
+  typedef AABB_primitive< IndexIterator,
+                          internal::Triangle_3_from_index_range_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>,
+                          internal::Point_from_indexed_triangle_3_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>,
+                          Tag_true,
+                          CacheDatum > Base;
+public:
+  ///constructor from an iterator
+  AABB_indexed_triangle_primitive_3(IndexIterator it, PointRange&) : Base(it) {}
+
+  /// \internal
+  static typename Base::Shared_data construct_shared_data(PointRange &range, PointMap pmap = PointMap()) {
+    return std::make_pair(
+      internal::Triangle_3_from_index_range_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>(range.begin(), pmap),
+      internal::Point_from_indexed_triangle_3_iterator_property_map<GeomTraits, IndexIterator, typename PointRange::iterator, PointMap>(range.begin(), pmap));
+  }
+};
+
+}  // end namespace CGAL
+
+#endif // CGAL_AABB_INDEXED_TRIANGLE_PRIMITIVE_3_H_
diff --git a/AABB_tree/include/CGAL/AABB_tree.h b/AABB_tree/include/CGAL/AABB_tree.h
index 6382bef9976f..3ef9601bdb46 100644
--- a/AABB_tree/include/CGAL/AABB_tree.h
+++ b/AABB_tree/include/CGAL/AABB_tree.h
@@ -32,6 +32,10 @@
 #include <CGAL/mutex.h>
 #endif
 
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/tbb.h>
+#endif
+
 /// \file AABB_tree.h
 
 namespace CGAL {
@@ -139,23 +143,35 @@ namespace CGAL {
     /// after one or more calls to `insert()`.
     /// This procedure is called implicitly at the first call to a query member function.
     /// An explicit call to `build()` must be made to ensure that the next call to
-    /// a query function will not trigger the construction of the data structure.
+    /// a query function will not trigger the construction of the data structure and/or
+    /// allow the tree to be built in parallel if supported.
+    /// `ConcurrencyTag` enables sequential versus parallel algorithm. Possible values are Sequential_tag, Parallel_tag, and Parallel_if_available_tag.
     /// A call to \link AABBTraits::set_shared_data `AABBTraits::set_shared_data(t...)`\endlink
-    // is made using the internally stored traits.
+    /// is made using the internally stored traits.
     /// This procedure has a complexity of \cgalBigO{n log(n)}, where \f$n\f$ is the number of
     /// primitives of the tree.
-    template<typename ... T>
+    template<typename ConcurrencyTag=Sequential_tag, typename ... T>
     void build(T&& ...);
 #ifndef DOXYGEN_RUNNING
+    template<typename ConcurrencyTag=Sequential_tag>
     void build();
 
     /// triggers the (re)construction of the tree similarly to a call to `build()`
     /// but the traits functors `Compute_bbox` and `Split_primitives` are ignored
     /// and `compute_bbox` and `split_primitives` are used instead.
-    template <class ComputeBbox, class SplitPrimitives>
+    template <class ConcurrencyTag=Sequential_tag, class ComputeBbox, class SplitPrimitives>
     void custom_build(const ComputeBbox& compute_bbox,
                       const SplitPrimitives& split_primitives);
 #endif
+    /// @private
+    template<typename ConcurrencyTag=Sequential_tag>
+    void partial_build(const size_t cutoff = 200000);
+
+    /// @private
+    template <class ConcurrencyTag=Sequential_tag, class ComputeBbox, class SplitPrimitives>
+    void custom_partial_build(const size_t cutoff,
+                              const ComputeBbox& compute_bbox,
+                              const SplitPrimitives& split_primitives);
     ///@}
 
     /// \name Operations
@@ -163,7 +179,7 @@ namespace CGAL {
 
     /// is equivalent to calling `clear()`, \link insert(InputIterator, InputIterator, T&&...) `insert(first,last,t...)`\endlink,
     // and `build()`
-    template<typename ConstPrimitiveIterator,typename ... T>
+    template<typename ConcurrencyTag=Sequential_tag, typename ConstPrimitiveIterator,typename ... T>
     void rebuild(ConstPrimitiveIterator first, ConstPrimitiveIterator beyond,T&& ...);
 
     /// adds a sequence of primitives to the set of primitives of the AABB tree.
@@ -228,8 +244,9 @@ namespace CGAL {
       set_primitive_data_impl(CGAL::Boolean_tag<internal::Has_nested_type_Shared_data<Primitive>::value>(),std::forward<T>(t)...);
     }
 
+    template<typename ConcurrencyTag=Sequential_tag>
     bool build_kd_tree();
-    template<typename ConstPointIterator>
+    template<typename ConcurrencyTag=Sequential_tag, typename ConstPointIterator>
     bool build_kd_tree(ConstPointIterator first, ConstPointIterator beyond);
 public:
 
@@ -439,6 +456,7 @@ namespace CGAL {
     /// constructs the internal search tree from
     /// a point set taken on the internal primitives
     /// returns `true` iff successful memory allocation
+    template<typename ConcurrencyTag=Sequential_tag>
     bool accelerate_distance_queries();
     /// turns off the usage of the internal search tree and clears it if it was already constructed.
     void do_not_accelerate_distance_queries();
@@ -451,11 +469,11 @@ namespace CGAL {
     /// is needed to update the search tree.
     /// \tparam ConstPointIterator is an iterator with
     /// value type `Point_and_primitive_id`.
-    template<typename ConstPointIterator>
+    template<typename ConcurrencyTag=Sequential_tag, typename ConstPointIterator>
     bool accelerate_distance_queries(ConstPointIterator first, ConstPointIterator beyond)
     {
       m_use_default_search_tree = false;
-      return build_kd_tree(first,beyond);
+      return build_kd_tree<ConcurrencyTag>(first,beyond);
     }
 
     /// returns the minimum squared distance between the query point
@@ -569,14 +587,36 @@ namespace CGAL {
      *
      * [first,beyond[ is the range of primitives to be added to the tree.
      */
-    template<typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
+    template<typename ConcurrencyTag=Sequential_tag, typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
     void expand(Node& node,
+                std::size_t node_index,
                 ConstPrimitiveIterator first,
                 ConstPrimitiveIterator beyond,
                 const std::size_t range,
                 const ComputeBbox& compute_bbox,
                 const SplitPrimitives& split_primitives);
 
+    /// @private
+    template<typename ConcurrencyTag=Sequential_tag, typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
+    void partial_expand(Node& node,
+                        std::size_t node_index,
+                        ConstPrimitiveIterator first,
+                        ConstPrimitiveIterator beyond,
+                        const std::size_t range,
+                        const std::size_t cutoff,
+                        const ComputeBbox& compute_bbox,
+                        const SplitPrimitives& split_primitives);
+
+  public:
+    /// @private
+    std::pair<typename std::vector<Primitive>::iterator, typename std::vector<Primitive>::iterator>
+      partial_node_to_primitives_iterator(const Node& node){
+      const Primitive* begin = std::addressof(node.left_data());
+      const Primitive* end = std::addressof(node.right_data());
+      return std::make_pair(m_primitives.begin() + (begin - m_primitives.data()),
+                            m_primitives.begin() + (end   - m_primitives.data()));
+    }
+
   public:
     // returns a point which must be on one primitive
     Point_and_primitive_id any_reference_point_and_id() const
@@ -629,9 +669,9 @@ namespace CGAL {
     Primitives m_primitives;
     // tree nodes. first node is the root node
     std::vector<Node> m_nodes;
-    #ifdef CGAL_HAS_THREADS
+#ifdef CGAL_HAS_THREADS
     mutable CGAL_MUTEX build_mutex; // mutex used to protect const calls inducing build() and build_kd_tree()
-    #endif
+#endif
   public:
     const Node* root_node() const {
       CGAL_assertion(size() > 1);
@@ -651,11 +691,15 @@ namespace CGAL {
       return std::addressof(m_nodes[0]);
     }
 
-    Node& new_node()
+#ifndef CGAL_NO_DEPRECATED_CODE
+    //  Now unused by build() function
+    /// @private
+    CGAL_DEPRECATED Node& new_node()
     {
       m_nodes.emplace_back();
       return m_nodes.back();
     }
+#endif // CGAL_NO_DEPRECATED_CODE
   private:
     const Primitive& singleton_data() const {
       CGAL_assertion(size() == 1);
@@ -748,7 +792,7 @@ namespace CGAL {
 
   // Clears tree and insert a set of primitives
   template<typename Tr>
-  template<typename ConstPrimitiveIterator, typename ... T>
+  template<typename ConcurrencyTag, typename ConstPrimitiveIterator, typename ... T>
   void AABB_tree<Tr>::rebuild(ConstPrimitiveIterator first,
                               ConstPrimitiveIterator beyond,
                               T&& ... t)
@@ -759,11 +803,11 @@ namespace CGAL {
     // inserts primitives
     insert(first, beyond,std::forward<T>(t)...);
 
-    build();
+    build<ConcurrencyTag>();
   }
 
   template<typename Tr>
-  template<typename ... T>
+  template<typename ConcurrencyTag, typename ... T>
   void AABB_tree<Tr>::build(T&& ... t)
   {
     set_shared_data(std::forward<T>(t)...);
@@ -783,16 +827,22 @@ namespace CGAL {
 #endif
   }
 
+
   template<typename Tr>
-  template<typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
+  template<typename ConcurrencyTag, typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
   void
   AABB_tree<Tr>::expand(Node& node,
+                        std::size_t node_index,
                         ConstPrimitiveIterator first,
                         ConstPrimitiveIterator beyond,
                         const std::size_t range,
                         const ComputeBbox& compute_bbox,
                         const SplitPrimitives& split_primitives)
   {
+    // TODO refined this hardcode value
+#ifdef CGAL_LINKED_WITH_TBB
+    const std::size_t cutoff_parallel_call = 30000; // min size for parallel call
+#endif
     node.set_bbox(compute_bbox(first, beyond));
 
     // sort primitives along longest axis aabb
@@ -804,29 +854,49 @@ namespace CGAL {
       node.set_children(*first, *(first+1));
       break;
     case 3:
-      node.set_children(*first, new_node());
-      expand(node.right_child(), first+1, beyond, 2, compute_bbox, split_primitives);
+      node.set_children(*first, m_nodes[node_index+1]);
+      expand(node.right_child(), node_index+1, first+1, beyond, 2, compute_bbox, split_primitives);
       break;
     default:
       const std::size_t new_range = range/2;
-      node.set_children(new_node(), new_node());
-      expand(node.left_child(), first, first + new_range, new_range, compute_bbox, split_primitives);
-      expand(node.right_child(), first + new_range, beyond, range - new_range, compute_bbox, split_primitives);
+      node.set_children(m_nodes[node_index+1], m_nodes[node_index+new_range]);
+#ifdef CGAL_LINKED_WITH_TBB
+      if constexpr(std::is_same_v<Parallel_tag, ConcurrencyTag>)
+      {
+        if(range > cutoff_parallel_call){
+          oneapi::tbb::task_group tg;
+          tg.run([&]{
+                  expand<ConcurrencyTag>(node.left_child(), node_index+1, first, first + new_range, new_range, compute_bbox, split_primitives); }
+                );
+          expand<ConcurrencyTag>(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, compute_bbox, split_primitives);
+          tg.wait();
+        } else {
+          expand(node.left_child(), node_index+1, first, first + new_range, new_range, compute_bbox, split_primitives);
+          expand(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, compute_bbox, split_primitives);
+        }
+      }
+      else
+#endif
+      {
+        expand(node.left_child(), node_index+1, first, first + new_range, new_range, compute_bbox, split_primitives);
+        expand(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, compute_bbox, split_primitives);
+      }
     }
   }
 
 
   // Build the data structure, after calls to insert(..)
   template<typename Tr>
+  template <class ConcurrencyTag>
   void AABB_tree<Tr>::build()
   {
-    custom_build(m_traits.compute_bbox_object(),
-                 m_traits.split_primitives_object());
+    custom_build<ConcurrencyTag>(m_traits.compute_bbox_object(),
+                                  m_traits.split_primitives_object());
   }
 #ifndef DOXYGEN_RUNNING
   // Build the data structure, after calls to insert(..)
   template<typename Tr>
-  template <class ComputeBbox, class SplitPrimitives>
+  template <class ConcurrencyTag, class ComputeBbox, class SplitPrimitives>
   void AABB_tree<Tr>::custom_build(
     const ComputeBbox& compute_bbox,
     const SplitPrimitives& split_primitives)
@@ -836,14 +906,15 @@ namespace CGAL {
     if(m_primitives.size() > 1) {
 
       // allocates tree nodes
-      m_nodes.reserve(m_primitives.size()-1);
+      m_nodes.resize(m_primitives.size()-1);
 
       // constructs the tree
-      expand(new_node(),
-             m_primitives.begin(), m_primitives.end(),
-             m_primitives.size(),
-             compute_bbox,
-             split_primitives);
+      expand<ConcurrencyTag>(m_nodes[0],
+                              0,
+                              m_primitives.begin(), m_primitives.end(),
+                              m_primitives.size(),
+                              compute_bbox,
+                              split_primitives);
     }
 #ifdef CGAL_HAS_THREADS
     m_atomic_need_build.store(false, std::memory_order_release); // in case build() is triggered by a call to root_node()
@@ -855,6 +926,7 @@ namespace CGAL {
   // constructs the search KD tree from given points
   // to accelerate the distance queries
   template<typename Tr>
+  template<typename ConcurrencyTag>
   bool AABB_tree<Tr>::build_kd_tree()
   {
     // iterate over primitives to get reference points on them
@@ -864,18 +936,20 @@ namespace CGAL {
       points.push_back( Point_and_primitive_id( Helper::get_reference_point(p, m_traits), p.id() ) );
 
     // clears current KD tree
-    return build_kd_tree(points.begin(), points.end());
+    return build_kd_tree<ConcurrencyTag>(points.begin(), points.end());
   }
 
   // constructs the search KD tree from given points
   // to accelerate the distance queries
   template<typename Tr>
-  template<typename ConstPointIterator>
+  template<typename ConcurrencyTag, typename ConstPointIterator>
   bool AABB_tree<Tr>::build_kd_tree(ConstPointIterator first,
                                     ConstPointIterator beyond)
   {
     clear_search_tree();
-    m_p_search_tree = std::make_unique<const Search_tree>(first, beyond);
+    std::unique_ptr<Search_tree> p_search_tree = std::make_unique<Search_tree>(first, beyond);
+    p_search_tree->template build<ConcurrencyTag>();
+    m_p_search_tree = std::move(p_search_tree);
 #ifdef CGAL_HAS_THREADS
       m_atomic_search_tree_constructed.store(true, std::memory_order_release); // in case build_kd_tree() is triggered by a call to best_hint()
 #else
@@ -884,6 +958,97 @@ namespace CGAL {
       return true;
   }
 
+  template<typename Tr>
+  template <class ConcurrencyTag>
+  void AABB_tree<Tr>::partial_build(const std::size_t cutoff)
+  {
+    custom_partial_build<ConcurrencyTag>(cutoff,
+                                        m_traits.compute_bbox_object(),
+                                        m_traits.split_primitives_object());
+  }
+
+  // Build the data structure, after calls to insert(..)
+  template<typename Tr>
+  template <class ConcurrencyTag, class ComputeBbox, class SplitPrimitives>
+  void AABB_tree<Tr>::custom_partial_build(
+    const std::size_t cutoff,
+    const ComputeBbox& compute_bbox,
+    const SplitPrimitives& split_primitives)
+  {
+    clear_nodes();
+
+    if(m_primitives.size() > 1) {
+
+      // allocates tree nodes
+      m_nodes.resize(m_primitives.size()-1);
+
+      // constructs the tree
+      partial_expand<ConcurrencyTag>(m_nodes[0],
+                                      0,
+                                      m_primitives.begin(), m_primitives.end(),
+                                      m_primitives.size(),
+                                      cutoff,
+                                      compute_bbox,
+                                      split_primitives);
+    }
+#ifdef CGAL_HAS_THREADS
+    m_atomic_need_build.store(false, std::memory_order_release); // in case build() is triggered by a call to root_node()
+#else
+    m_need_build = false;
+#endif
+  }
+
+  template<typename Tr>
+  template<typename ConcurrencyTag, typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
+  void
+  AABB_tree<Tr>::partial_expand(Node& node,
+                                std::size_t node_index,
+                                ConstPrimitiveIterator first,
+                                ConstPrimitiveIterator beyond,
+                                const std::size_t range,
+                                const std::size_t cutoff,
+                                const ComputeBbox& compute_bbox,
+                                const SplitPrimitives& split_primitives)
+  {
+#ifdef CGAL_LINKED_WITH_TBB
+    const std::size_t cutoff_parallel_call = 30000; // min size for parallel call
+#endif
+    node.set_bbox(compute_bbox(first, beyond));
+
+    if(range < cutoff)
+    {
+      node.set_children(*first, *beyond);
+    }
+    else
+    {
+      // sort primitives along longest axis aabb
+      split_primitives(first, beyond, node.bbox());
+      const std::size_t new_range = range/2;
+      node.set_children(m_nodes[node_index+1], m_nodes[node_index+new_range]);
+#ifdef CGAL_LINKED_WITH_TBB
+      if constexpr(std::is_same_v<Parallel_tag, ConcurrencyTag>)
+      {
+        if(range > cutoff_parallel_call){
+          oneapi::tbb::task_group tg;
+          tg.run([&]{
+                  partial_expand<ConcurrencyTag>(node.left_child(), node_index+1, first, first + new_range, new_range, cutoff, compute_bbox, split_primitives); }
+                );
+          partial_expand<ConcurrencyTag>(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, cutoff, compute_bbox, split_primitives);
+          tg.wait();
+        } else {
+          partial_expand(node.left_child(), node_index+1, first, first + new_range, new_range, cutoff, compute_bbox, split_primitives);
+          partial_expand(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, cutoff, compute_bbox, split_primitives);
+        }
+      }
+      else
+#endif
+      {
+        partial_expand(node.left_child(), node_index+1, first, first + new_range, new_range, cutoff, compute_bbox, split_primitives);
+        partial_expand(node.right_child(), node_index+new_range, first + new_range, beyond, range - new_range, cutoff, compute_bbox, split_primitives);
+      }
+    }
+  }
+
   template<typename Tr>
   void AABB_tree<Tr>::do_not_accelerate_distance_queries()
   {
@@ -893,11 +1058,12 @@ namespace CGAL {
 
   // constructs the search KD tree from internal primitives
   template<typename Tr>
+  template<typename ConcurrencyTag>
   bool AABB_tree<Tr>::accelerate_distance_queries()
   {
     m_use_default_search_tree = true;
     if(m_primitives.empty()) return true;
-    return build_kd_tree();
+    return build_kd_tree<ConcurrencyTag>();
   }
 
   template<typename Tr>
diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_search_tree.h b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_search_tree.h
index f44130eda28b..9bc7f314d285 100644
--- a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_search_tree.h
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_search_tree.h
@@ -58,7 +58,12 @@ struct AABB_search_tree
       ++begin;
     }
     m_tree.insert(points.begin(), points.end());
-    m_tree.build();
+  }
+
+  template <class ConcurrencyTag=Sequential_tag>
+  void build()
+  {
+    m_tree.template build<ConcurrencyTag>();
   }
 
   template <typename Point>
diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_traversal_traits.h b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_traversal_traits.h
index 11fe903df932..fd3d6b694196 100644
--- a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_traversal_traits.h
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_traversal_traits.h
@@ -15,7 +15,7 @@
 
 #include <CGAL/license/AABB_tree.h>
 
-
+#include <CGAL/AABB_tree/internal/Primitive_helper.h>
 #include <CGAL/AABB_tree/internal/AABB_node.h>
 #include <optional>
 
@@ -183,6 +183,96 @@ class Listing_primitive_traits
   const AABBTraits& m_traits;
 };
 
+/**
+ * @class Listing_primitive_traits_with_transformation
+ */
+template<typename AABBTraits, typename Kernel, typename Query, typename Output_iterator>
+class Listing_primitive_traits_with_transformation
+{
+  typedef typename AABBTraits::FT FT;
+  typedef typename AABBTraits::Point Point;
+  typedef typename AABBTraits::Primitive Primitive;
+  typedef typename AABBTraits::Bounding_box Bounding_box;
+  typedef typename AABBTraits::Primitive::Id Primitive_id;
+  typedef typename AABBTraits::Point_and_primitive_id Point_and_primitive_id;
+  typedef typename AABBTraits::Object_and_primitive_id Object_and_primitive_id;
+  typedef ::CGAL::AABB_node<AABBTraits> Node;
+
+  Bbox_3 compute_transformed_bbox(const Bbox_3& bbox) const
+  {
+    return internal::compute_transformed_bbox(m_transfo, bbox, m_has_rotation);
+  }
+
+public:
+  Listing_primitive_traits_with_transformation(Output_iterator out_it,const AABBTraits& traits, const Aff_transformation_repC3<Kernel> &transfo)
+    : m_out_it(out_it), m_traits(traits), m_transfo(transfo), m_has_rotation(has_rotation(transfo))
+  {}
+
+  bool go_further() const { return true; }
+
+  void intersection(const Query& query, const Primitive& primitive)
+  {
+    auto datum = internal::Primitive_helper<AABBTraits>::get_datum(primitive, m_traits);
+    if( CGAL::do_intersect(query, datum.transform(Aff_transformation_3<Kernel>(m_transfo))) )
+      *m_out_it++ = primitive.id();
+  }
+
+  bool do_intersect(const Query& query, const Node& node) const
+  {
+#if 1
+    return m_traits.do_intersect_object()(query, compute_transformed_bbox(node.bbox()));
+#else
+    return Convex_hull_3::do_intersect(query, node.bbox(), parameters::transformation(Aff_transformation_3<Kernel>(m_transfo)))
+#endif
+  }
+
+private:
+  Output_iterator m_out_it;
+  const AABBTraits m_traits;
+  const Aff_transformation_repC3<Kernel>& m_transfo;
+  bool m_has_rotation;
+};
+
+/**
+ * @class Listing_distinct_primitive_traits
+ * used by `all_pairs_of_intersecting_primitives()` to avoid report `(i, i)` and twice `(i, j)`.
+ */
+template<typename AABBTraits, typename Output_iterator>
+class Listing_distinct_primitive_traits
+{
+  typedef typename AABBTraits::FT FT;
+  typedef typename AABBTraits::Point Point;
+  typedef typename AABBTraits::Primitive Primitive;
+  typedef typename AABBTraits::Bounding_box Bounding_box;
+  typedef typename AABBTraits::Primitive::Id Primitive_id;
+  typedef typename AABBTraits::Point_and_primitive_id Point_and_primitive_id;
+  typedef typename AABBTraits::Object_and_primitive_id Object_and_primitive_id;
+  typedef ::CGAL::AABB_node<AABBTraits> Node;
+
+public:
+  Listing_distinct_primitive_traits(Output_iterator out_it, const AABBTraits& traits)
+    : m_out_it(out_it), m_traits(traits) {}
+
+  constexpr bool go_further() const { return true; }
+
+  void intersection(const Primitive& query, const Primitive& primitive)
+  {
+    if( query.id()<primitive.id() && m_traits.do_intersect_object()(internal::Primitive_helper<AABBTraits>::get_datum(query, m_traits), primitive) )
+    {
+      *m_out_it++ = primitive.id();
+    }
+  }
+
+  bool do_intersect(const Primitive& query, const Node& node) const
+  {
+    return m_traits.do_intersect_object()(query, node.bbox());
+  }
+
+private:
+  Output_iterator m_out_it;
+  const AABBTraits& m_traits;
+};
+
 
 /**
  * @class First_primitive_traits
@@ -270,6 +360,57 @@ class Do_intersect_traits
   const AABBTraits& m_traits;
 };
 
+/**
+ * @class Do_intersect_traits_with_transformation
+ */
+template<typename AABBTraits, typename Kernel, typename Query>
+class Do_intersect_traits_with_transformation
+{
+  typedef typename AABBTraits::FT FT;
+  typedef typename AABBTraits::Point Point;
+  typedef typename AABBTraits::Primitive Primitive;
+  typedef typename AABBTraits::Bounding_box Bounding_box;
+  typedef typename AABBTraits::Primitive::Id Primitive_id;
+  typedef typename AABBTraits::Point_and_primitive_id Point_and_primitive_id;
+  typedef typename AABBTraits::Object_and_primitive_id Object_and_primitive_id;
+  typedef ::CGAL::AABB_node<AABBTraits> Node;
+
+  Bbox_3 compute_transformed_bbox(const Bbox_3& bbox) const
+  {
+    return compute_transformed_bbox(m_transfo, bbox, m_has_rotation);
+  }
+
+public:
+  Do_intersect_traits_with_transformation(const AABBTraits& traits, const Aff_transformation_3<Kernel> &transfo)
+    : m_is_found(false), m_traits(traits), m_transfo(transfo), m_has_rotation(has_rotation(transfo))
+  {}
+
+  bool go_further() const { return !m_is_found; }
+
+  void intersection(const Query& query, const Primitive& primitive)
+  {
+    if( m_traits.do_intersect_object()(query, internal::Primitive_helper<AABBTraits>::get_datum(primitive, m_traits).transform(Aff_transformation_3<Kernel>(m_transfo))) )
+      m_is_found = true;
+  }
+
+  bool do_intersect(const Query& query, const Node& node) const
+  {
+#if 1
+    return m_traits.do_intersect_object()(query, compute_transformed_bbox(node.bbox()));
+#else
+    return Convex_hull_3::do_intersect(query, node.bbox(), parameters::transformation(m_transfo))
+#endif
+  }
+
+  bool is_intersection_found() const { return m_is_found; }
+
+private:
+  bool m_is_found;
+  const AABBTraits m_traits;
+  Aff_transformation_3<Kernel> m_transfo;
+  bool m_has_rotation;
+};
+
 
 /**
  * @class Projection_traits
diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal.h b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal.h
new file mode 100644
index 000000000000..1c8c9e8a683a
--- /dev/null
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal.h
@@ -0,0 +1,197 @@
+// Copyright (c) 2026  Geometry Factory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s) : Léo Valque
+
+#ifndef CGAL_AABB_TWO_TREES_TRAVERSAL_H
+#define CGAL_AABB_TWO_TREES_TRAVERSAL_H
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/AABB_tree/internal/AABB_node.h>
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+
+namespace CGAL {
+
+namespace internal { namespace AABB_tree {
+
+template <bool in_order = true,
+          typename ConcurrencyTag = Sequential_tag,
+          typename AABBTraits_A,
+          typename AABBTraits_B,
+          typename TwoTreeTraversalTraits>
+void two_trees_traversal(const ::CGAL::AABB_node<AABBTraits_A>& node_A,
+                        const ::CGAL::AABB_node<AABBTraits_B>& node_B,
+                        const std::size_t nb_primitives_A,
+                        const std::size_t nb_primitives_B,
+                        TwoTreeTraversalTraits& traversal_traits)
+{
+#if CGAL_LINKED_WITH_TBB
+  const std::size_t cutoff_parallel_call = 100000;
+#endif
+  auto recursive_call = [](const auto &node_A, const auto &node_B, const std::size_t nb_primitives_A, const std::size_t nb_primitives_B, auto &traversal_traits){
+    if(traversal_traits.prefer_A_for_next_step(node_A, node_B, nb_primitives_A, nb_primitives_B))
+      two_trees_traversal< in_order, ConcurrencyTag>(node_A, node_B, nb_primitives_A, nb_primitives_B, traversal_traits);
+    else
+      two_trees_traversal<!in_order, ConcurrencyTag>(node_B, node_A, nb_primitives_B, nb_primitives_A, traversal_traits);
+  };
+  switch(nb_primitives_A)
+  {
+  case 2:
+  {
+    if constexpr(in_order){
+      traversal_traits.intersection(node_A.left_data(), node_B, nb_primitives_B);
+      traversal_traits.intersection(node_A.right_data(), node_B, nb_primitives_B);
+    } else {
+      traversal_traits.intersection(node_B, nb_primitives_B, node_A.left_data());
+      traversal_traits.intersection(node_B, nb_primitives_B, node_A.right_data());
+    }
+    break;
+  }
+  case 3:
+  {
+    if constexpr(in_order)
+      traversal_traits.intersection(node_A.left_data(), node_B, nb_primitives_B);
+    else
+      traversal_traits.intersection(node_B, nb_primitives_B, node_A.left_data());
+
+    if( traversal_traits.do_intersect(node_A.right_child(), node_B) )
+      two_trees_traversal<!in_order>(node_B, node_A.right_child(), nb_primitives_B, 2, traversal_traits);
+    break;
+  }
+  default:
+  {
+    bool do_intersect_left  = traversal_traits.do_intersect(node_A.left_child(), node_B);
+    bool do_intersect_right = traversal_traits.do_intersect(node_A.right_child(), node_B);
+#if CGAL_LINKED_WITH_TBB
+    if constexpr(std::is_same_v<ConcurrencyTag, Parallel_tag>)
+    {
+      if(do_intersect_left && do_intersect_right && nb_primitives_A > cutoff_parallel_call && nb_primitives_B > cutoff_parallel_call)
+      {
+        oneapi::tbb::task_group tg;
+        tg.run([&]{
+                recursive_call(node_A.left_child(), node_B, nb_primitives_A/2, nb_primitives_B, traversal_traits);}
+              );
+        recursive_call(node_A.right_child(), node_B, nb_primitives_A - nb_primitives_A/2, nb_primitives_B, traversal_traits);
+        tg.wait();
+      }
+      else
+      {
+        if( do_intersect_left )
+          recursive_call(node_A.left_child(), node_B, nb_primitives_A/2, nb_primitives_B, traversal_traits);
+        if( traversal_traits.go_further() && do_intersect_right )
+          recursive_call(node_A.right_child(), node_B, nb_primitives_A - nb_primitives_A/2, nb_primitives_B, traversal_traits);
+      }
+    }
+    else
+#endif
+    {
+      if( do_intersect_left )
+        recursive_call(node_A.left_child(), node_B, nb_primitives_A/2, nb_primitives_B, traversal_traits);
+      if( traversal_traits.go_further() && do_intersect_right )
+        recursive_call(node_A.right_child(), node_B, nb_primitives_A - nb_primitives_A/2, nb_primitives_B, traversal_traits);
+    }
+  }} // switch end
+}
+
+template<typename ConcurrencyTag = Sequential_tag,
+         typename Tree_A,
+         typename Tree_B,
+         typename TwoTreeTraversalTraits>
+void two_trees_traversal(const Tree_A& tree_A,
+                        const Tree_B& tree_B,
+                        TwoTreeTraversalTraits &traits)
+{
+  CGAL_precondition(tree_A.size() != 0 && tree_B.size() != 0);
+  two_trees_traversal<true, ConcurrencyTag>(*tree_A.root_node(), *tree_B.root_node(), tree_A.size(), tree_B.size(), traits);
+}
+
+namespace experimental{
+
+template <bool in_order = true,
+          typename ConcurrencyTag = Sequential_tag,
+          typename AABBTraits_A,
+          typename AABBTraits_B,
+          typename TwoTreeTraversalTraits>
+void two_trees_partial_traversal(const ::CGAL::AABB_node<AABBTraits_A>& node_A,
+                                const ::CGAL::AABB_node<AABBTraits_B>& node_B,
+                                const std::size_t nb_primitives_A,
+                                const std::size_t nb_primitives_B,
+                                const std::size_t cutoff,
+                                TwoTreeTraversalTraits& traversal_traits)
+{
+#if CGAL_LINKED_WITH_TBB
+  const std::size_t cutoff_parallel_call = 100000;
+#endif
+  if(nb_primitives_A < cutoff && nb_primitives_B < cutoff)
+  {
+    if constexpr(in_order)
+      traversal_traits.intersection(node_A, node_B);
+    else
+      traversal_traits.intersection(node_B, node_A);
+  }
+  else if(nb_primitives_A < cutoff && nb_primitives_B < cutoff)
+  {
+    two_trees_partial_traversal<!in_order>(node_B, node_A, nb_primitives_B, nb_primitives_A, cutoff, traversal_traits);
+  }
+  else
+  {
+    // TODO current strategy is we swap at each step. Try largest tree first or largest bbox first
+    bool do_intersect_left  = traversal_traits.do_intersect(node_A.left_child(), node_B);
+    bool do_intersect_right = traversal_traits.do_intersect(node_A.right_child(), node_B);
+#if CGAL_LINKED_WITH_TBB
+    if constexpr(std::is_same_v<ConcurrencyTag, Parallel_tag>)
+    {
+      if(do_intersect_left && do_intersect_right && nb_primitives_A > cutoff_parallel_call && nb_primitives_B > cutoff_parallel_call)
+      {
+        oneapi::tbb::task_group tg;
+        tg.run([&]{
+                two_trees_partial_traversal<!in_order, ConcurrencyTag>(node_B, node_A.left_child(), nb_primitives_B, nb_primitives_A/2, cutoff, traversal_traits);}
+              );
+        two_trees_partial_traversal<!in_order, ConcurrencyTag>(node_B, node_A.right_child(), nb_primitives_B, nb_primitives_A-nb_primitives_A/2, cutoff, traversal_traits);
+        tg.wait();
+      }
+      else
+      {
+        if( do_intersect_left )
+          two_trees_partial_traversal<!in_order>(node_B, node_A.left_child(), nb_primitives_B, nb_primitives_A/2, cutoff, traversal_traits);
+        if( traversal_traits.go_further() && do_intersect_right )
+          two_trees_partial_traversal<!in_order>(node_B, node_A.right_child(), nb_primitives_B, nb_primitives_A-nb_primitives_A/2, cutoff, traversal_traits);
+      }
+    }
+    else
+#endif
+    {
+      if( do_intersect_left )
+        two_trees_partial_traversal<!in_order>(node_B, node_A.left_child(), nb_primitives_B, nb_primitives_A/2, cutoff, traversal_traits);
+      if( traversal_traits.go_further() && do_intersect_right )
+        two_trees_partial_traversal<!in_order>(node_B, node_A.right_child(), nb_primitives_B,  nb_primitives_A-nb_primitives_A/2, cutoff, traversal_traits);
+    }
+  }
+}
+
+template<typename ConcurrencyTag = Sequential_tag,
+         typename Tree_A,
+         typename Tree_B,
+         typename TwoTreeTraversalTraits>
+void two_trees_partial_traversal(const Tree_A& tree_A,
+                                const Tree_B& tree_B,
+                                const std::size_t cutoff,
+                                TwoTreeTraversalTraits &traits)
+{
+  CGAL_precondition(tree_A.size() != 0 && tree_B.size() != 0);
+  two_trees_partial_traversal<true, ConcurrencyTag>(*tree_A.root_node(), *tree_B.root_node(), tree_A.size(), tree_B.size(), cutoff, traits);
+}
+
+} // end of namespace experimental
+
+}}} // end of namespace CGAL::internal::AABB_tree
+
+#endif // CGAL_AABB_TRAVERSAL_TRAITS_H
diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal_traits.h b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal_traits.h
new file mode 100644
index 000000000000..159f30718ce0
--- /dev/null
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/AABB_two_trees_traversal_traits.h
@@ -0,0 +1,376 @@
+// Copyright (c) 2026  Geometry Factory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s) : Léo Valque
+
+#ifndef CGAL_AABB_TWO_TREES_TRAVERSAL_TRAITS_H
+#define CGAL_AABB_TWO_TREES_TRAVERSAL_TRAITS_H
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/AABB_tree/internal/Primitive_helper.h>
+#include <CGAL/AABB_tree/internal/AABB_node.h>
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+#include <optional>
+
+namespace CGAL {
+
+namespace internal { namespace AABB_tree {
+
+template<bool in_order, typename Value, typename OutputIterator>
+class Wrap_output_iterator
+{
+  Value first;
+  OutputIterator out;
+public:
+  Wrap_output_iterator(Value first_, OutputIterator out_): first(first_), out(out_){}
+
+  Wrap_output_iterator& operator=(Value second){
+    if constexpr(in_order)
+      out = std::make_pair(first, second);
+    else
+      out = std::make_pair(second, first);
+    return *this;
+  }
+  Wrap_output_iterator& operator*(){ return *this; }
+  Wrap_output_iterator& operator++(){ ++out; return *this; }
+  Wrap_output_iterator& operator++(int){ /*auto tmp = *this;*/ ++out; return *this; }
+  Wrap_output_iterator& operator+(int d){ out += d; return *this; }
+};
+
+template<typename AABBTraits1, typename AABBTraits2, typename OutputIterator>
+class Two_trees_listing_intersecting_primitives_traits
+{
+  typedef typename AABBTraits1::Primitive Primitive1;
+  typedef typename AABBTraits2::Primitive Primitive2;
+  typedef ::CGAL::AABB_node<AABBTraits1> Node1;
+  typedef ::CGAL::AABB_node<AABBTraits2> Node2;
+
+public:
+  Two_trees_listing_intersecting_primitives_traits(const AABBTraits1& traits1, const AABBTraits2& traits2, OutputIterator out_)
+    : m_traits1(traits1), m_traits2(traits2), out(out_)
+  {}
+
+  bool go_further() const {
+    return true;
+  }
+
+  // If true, the next step of traversal continue on A, if false, the next step of traversal is on B
+  template<typename Node_A, typename Node_B>
+  bool prefer_A_for_next_step(const Node_A& node_a, const Node_B& node_b, const std::size_t&, const std::size_t&) const {
+  #if 1
+    return node_a.bbox().squared_diagonal_length() > node_b.bbox().squared_diagonal_length();
+  #else
+    return false;
+  #endif
+  }
+
+  void intersection(const Primitive1& primitive1, const Node2& node2, std::size_t nb_primitives_2)
+  {
+    using Wrap_iterator = Wrap_output_iterator<true, typename Primitive1::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive1.id(), out);
+    Listing_primitive_traits<AABBTraits2, typename AABBTraits1::Primitive::Datum, Wrap_iterator> traits(wrap_out, m_traits2);
+    node2.traversal( internal::Primitive_helper<AABBTraits1>::get_datum(primitive1, m_traits1), traits, nb_primitives_2);
+  }
+
+  void intersection(const Node1& node1, std::size_t nb_primitives_1, const Primitive2& primitive2)
+  {
+    using Wrap_iterator= Wrap_output_iterator<false, typename Primitive2::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive2.id(), out);
+    Listing_primitive_traits<AABBTraits1, typename AABBTraits2::Primitive::Datum, Wrap_iterator> traits(wrap_out, m_traits1);
+    node1.traversal( internal::Primitive_helper<AABBTraits2>::get_datum(primitive2, m_traits2), traits, nb_primitives_1);
+  }
+
+  bool do_intersect(const Node1& node1, const Node2& node2) const
+  {
+    return do_overlap(node1.bbox(), node2.bbox());
+  }
+
+private:
+  const AABBTraits1& m_traits1;
+  const AABBTraits2& m_traits2;
+  OutputIterator out;
+};
+
+template<typename AABBTraits1, typename AABBTraits2, typename OutputIterator, typename Kernel>
+class Two_trees_listing_intersecting_primitives_traits_with_transformation
+{
+  typedef typename AABBTraits1::Primitive Primitive1;
+  typedef typename AABBTraits2::Primitive Primitive2;
+  typedef ::CGAL::AABB_node<AABBTraits1> Node1;
+  typedef ::CGAL::AABB_node<AABBTraits2> Node2;
+
+public:
+  Two_trees_listing_intersecting_primitives_traits_with_transformation(const AABBTraits1& traits1, const AABBTraits2& traits2,
+                                                                       OutputIterator out_,
+                                                                       const Aff_transformation_repC3<Kernel> &tr1, const Aff_transformation_repC3<Kernel> &tr2)
+    : m_traits1(traits1), m_traits2(traits2), out(out_), m_tr1(tr1), m_tr2(tr2), m_tr1_has_rotation(has_rotation(tr1)), m_tr2_has_rotation(has_rotation(tr2))
+  {}
+
+  bool go_further() const {
+    return true;
+  }
+
+  // If true, the next step of traversal continue on A, if false, the next step of traversal is on B
+  template<typename Node_A, typename Node_B>
+  bool prefer_A_for_next_step(const Node_A& node_a, const Node_B& node_b, const std::size_t&, const std::size_t&) const {
+  #if 1
+    return node_a.bbox().squared_diagonal_length() > node_b.bbox().squared_diagonal_length();
+  #else
+    return false;
+  #endif
+  }
+
+  void intersection(const Primitive1& primitive1, const Node2& node2, std::size_t nb_primitives_2)
+  {
+    using Wrap_iterator = Wrap_output_iterator<true, typename Primitive1::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive1.id(), out);
+    Listing_primitive_traits_with_transformation<AABBTraits2, Kernel, typename AABBTraits1::Primitive::Datum, Wrap_iterator> traits(wrap_out, m_traits2, m_tr2);
+    auto datum = internal::Primitive_helper<AABBTraits1>::get_datum(primitive1, m_traits1).transform(Aff_transformation_3<Kernel>(m_tr1));
+    node2.traversal(datum, traits, nb_primitives_2);
+  }
+
+  void intersection(const Node1& node1, std::size_t nb_primitives_1, const Primitive2& primitive2)
+  {
+    using Wrap_iterator= Wrap_output_iterator<false, typename Primitive2::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive2.id(), out);
+    Listing_primitive_traits_with_transformation<AABBTraits1, Kernel, typename AABBTraits2::Primitive::Datum, Wrap_iterator> traits(wrap_out, m_traits1, m_tr1);
+    auto datum = internal::Primitive_helper<AABBTraits1>::get_datum(primitive2, m_traits2).transform(Aff_transformation_3<Kernel>(m_tr2));
+    node1.traversal(datum, traits, nb_primitives_1);
+  }
+
+  bool do_intersect(const Node1& node1, const Node2& node2) const
+  {
+  // TODO Both option are slow
+  #if 1
+    return do_overlap(compute_transformed_bbox(m_tr1, node1.bbox(), m_tr1_has_rotation), compute_transformed_bbox(m_tr2, node2.bbox(), m_tr2_has_rotation));
+  #else
+    return Convex_hull_3::do_intersect(node1.bbox(), node2.bbox(), parameters::transformation(m_tr1), parameters::transformation(m_tr2));
+  #endif
+  }
+
+private:
+  const AABBTraits1& m_traits1;
+  const AABBTraits2& m_traits2;
+  OutputIterator out;
+  const Aff_transformation_repC3<Kernel>& m_tr1;
+  const Aff_transformation_repC3<Kernel>& m_tr2;
+  bool m_tr1_has_rotation, m_tr2_has_rotation;
+};
+
+template<typename AABBTraits, typename OutputIterator>
+class Listing_self_intersecting_primitives_traits
+{
+  typedef typename AABBTraits::Primitive Primitive;
+  typedef ::CGAL::AABB_node<AABBTraits> Node;
+
+public:
+  Listing_self_intersecting_primitives_traits(const AABBTraits& traits, OutputIterator out_)
+    : m_traits(traits), out(out_)
+  {}
+
+  bool go_further() const {
+    return true;
+  }
+
+  // If true, the next step of traversal continue on A, if false, the next step of traversal is on B
+  template<typename Node_A, typename Node_B>
+  bool prefer_A_for_next_step(const Node_A& node_a, const Node_B& node_b, const std::size_t&, const std::size_t&) const {
+    return node_a.bbox().squared_diagonal_length() > node_b.bbox().squared_diagonal_length();
+  }
+
+  void intersection(const Primitive& primitive1, const Node& node2, std::size_t nb_primitives_2)
+  {
+    // TODO Since we are in a symmetric case, maybe we can ignore this call
+    using Wrap_iterator = Wrap_output_iterator<true, typename Primitive::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive1.id(), out);
+    Listing_distinct_primitive_traits<AABBTraits, Wrap_iterator> traits(wrap_out, m_traits);
+    node2.traversal( primitive1, traits, nb_primitives_2);
+  }
+
+  void intersection(const Node& node1, std::size_t nb_primitives_1, const Primitive& primitive2)
+  {
+    using Wrap_iterator= Wrap_output_iterator<false, typename Primitive::Id, OutputIterator>;
+    Wrap_iterator wrap_out(primitive2.id(), out);
+    Listing_distinct_primitive_traits<AABBTraits, Wrap_iterator> traits(wrap_out, m_traits);
+    node1.traversal( primitive2, traits, nb_primitives_1);
+  }
+
+  bool do_intersect(const Node& node1, const Node& node2) const
+  {
+    return do_overlap(node1.bbox(), node2.bbox());
+  }
+
+private:
+  const AABBTraits& m_traits;
+  OutputIterator out;
+};
+
+template<typename AABBTraits1, typename AABBTraits2>
+class Two_trees_do_intersect_traits
+{
+  typedef typename AABBTraits1::Primitive Primitive1;
+  typedef typename AABBTraits2::Primitive Primitive2;
+  typedef ::CGAL::AABB_node<AABBTraits1> Node1;
+  typedef ::CGAL::AABB_node<AABBTraits2> Node2;
+
+public:
+  Two_trees_do_intersect_traits(const AABBTraits1& traits1, const AABBTraits2& traits2)
+    : m_traits1(traits1), m_traits2(traits2), m_is_found(false)
+  {}
+
+  bool go_further() const {
+    return !m_is_found;
+  }
+
+  // If true, the next step of traversal continues on A, if false, the next step of traversal is on B
+  template<typename Node_A, typename Node_B>
+  bool prefer_A_for_next_step(const Node_A& node_a, const Node_B& node_b, const std::size_t&, const std::size_t&) const {
+  #if 1
+    return node_a.bbox().squared_diagonal_length() > node_b.bbox().squared_diagonal_length();
+  #else
+    return false;
+  #endif
+  }
+
+  void intersection(const Primitive1& primitive1, const Node2& node2, std::size_t nb_primitives_2)
+  {
+    Do_intersect_traits<AABBTraits2, typename AABBTraits1::Primitive::Datum> traits(m_traits2);
+    node2.traversal( internal::Primitive_helper<AABBTraits1>::get_datum(primitive1, m_traits1), traits, nb_primitives_2);
+    if(traits.is_intersection_found())
+      m_is_found = true;
+  }
+
+  void intersection(const Node1& node1, std::size_t nb_primitives_1, const Primitive2& primitive2)
+  {
+    Do_intersect_traits<AABBTraits1, typename AABBTraits2::Primitive::Datum> traits(m_traits1);
+    node1.traversal( internal::Primitive_helper<AABBTraits2>::get_datum(primitive2, m_traits2), traits, nb_primitives_1);
+    if(traits.is_intersection_found())
+      m_is_found = true;
+  }
+
+  bool do_intersect(const Node1& node1, const Node2& node2) const
+  {
+    return do_overlap(node1.bbox(), node2.bbox());
+  }
+
+  bool is_intersection_found() const { return m_is_found; }
+
+private:
+  const AABBTraits1& m_traits1;
+  const AABBTraits2& m_traits2;
+  bool m_is_found;
+};
+
+template<typename AABBTraits1, typename AABBTraits2, typename Kernel>
+class Two_trees_do_intersect_traits_with_transformation
+{
+  typedef typename AABBTraits1::Primitive Primitive1;
+  typedef typename AABBTraits2::Primitive Primitive2;
+  typedef ::CGAL::AABB_node<AABBTraits1> Node1;
+  typedef ::CGAL::AABB_node<AABBTraits2> Node2;
+
+public:
+  Two_trees_do_intersect_traits_with_transformation(const AABBTraits1& traits1, const AABBTraits2& traits2, const Aff_transformation_repC3<Kernel> &tr1, const Aff_transformation_repC3<Kernel> &tr2)
+    : m_traits1(traits1), m_traits2(traits2), m_tr1(tr1), m_tr2(tr2), m_is_found(false)
+  {}
+
+  bool go_further() const {
+    return !m_is_found;
+  }
+
+  // If true, the next step of traversal continues on A, if false, the next step of traversal is on B
+  template<typename Node_A, typename Node_B>
+  bool prefer_A_for_next_step(const Node_A& node_a, const Node_B& node_b, const std::size_t&, const std::size_t&) const {
+  #if 1
+    return node_a.bbox().squared_diagonal_length() > node_b.bbox().squared_diagonal_length();
+  #else
+    return false;
+  #endif
+  }
+
+  void intersection(const Primitive1& primitive1, const Node2& node2, std::size_t nb_primitives_2)
+  {
+    Do_intersect_traits_with_transformation<AABBTraits2, typename AABBTraits1::Primitive::Datum, Kernel> traits(m_traits2, m_tr2);
+    auto datum = internal::Primitive_helper<AABBTraits1>::get_datum(primitive1, m_traits1).transform(m_tr1);
+    node2.traversal( datum, traits, nb_primitives_2);
+    if(traits.is_intersection_found())
+      m_is_found = true;
+  }
+
+  void intersection(const Node1& node1, std::size_t nb_primitives_1, const Primitive2& primitive2)
+  {
+    Do_intersect_traits_with_transformation<AABBTraits1, typename AABBTraits2::Primitive::Datum, Kernel> traits(m_traits1, m_tr1);
+    auto datum = internal::Primitive_helper<AABBTraits1>::get_datum(primitive2, m_traits2).transform(m_tr2);
+    node1.traversal( datum, traits, nb_primitives_1);
+    if(traits.is_intersection_found())
+      m_is_found = true;
+  }
+
+  bool do_intersect(const Node1& node1, const Node2& node2) const
+  {
+  #if 1
+    return do_overlap(compute_transformed_bbox(node1.bbox(), m_tr1), compute_transformed_bbox(node2.bbox(), m_tr2));
+  #else
+    return Convex_hull_3::do_intersect(node1.bbox(), node2.bbox(), parameters::transformation(m_tr1), parameters::transformation(m_tr2));
+  #endif
+  }
+
+  bool is_intersection_found() const { return m_is_found; }
+
+private:
+  const AABBTraits1& m_traits1;
+  const AABBTraits2& m_traits2;
+  const Aff_transformation_repC3<Kernel>& m_tr1;
+  const Aff_transformation_repC3<Kernel>& m_tr2;
+  bool m_is_found;
+};
+
+namespace experimental{
+
+template<typename AABBTraits1, typename AABBTraits2, typename OutputIterator>
+class Two_trees_intersecting_nodes_traits
+{
+  typedef typename AABBTraits1::Primitive Primitive1;
+  typedef typename AABBTraits2::Primitive Primitive2;
+  typedef ::CGAL::AABB_node<AABBTraits1> Node1;
+  typedef ::CGAL::AABB_node<AABBTraits2> Node2;
+
+public:
+  Two_trees_intersecting_nodes_traits(const AABBTraits1& traits1, const AABBTraits2& traits2, OutputIterator out_)
+    : m_traits1(traits1), m_traits2(traits2), out(out_)
+  {}
+
+  bool go_further() const {
+    return true;
+  }
+
+  void intersection(const Node1& node1, const Node2& node2)
+  {
+    *out++ = std::make_pair(&node1, &node2);
+  }
+
+  bool do_intersect(const Node1& node1, const Node2& node2) const
+  {
+    return do_overlap(node1.bbox(), node2.bbox());
+  }
+
+private:
+  const AABBTraits1& m_traits1;
+  const AABBTraits2& m_traits2;
+  OutputIterator out;
+};
+
+}
+
+
+}}} // end namespace CGAL::internal::AABB_tree
+
+#endif // CGAL_AABB_TRAVERSAL_TRAITS_H
diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/Primitive_helper.h b/AABB_tree/include/CGAL/AABB_tree/internal/Primitive_helper.h
index b8bfa8a199c8..cd889337c1b4 100644
--- a/AABB_tree/include/CGAL/AABB_tree/internal/Primitive_helper.h
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/Primitive_helper.h
@@ -14,7 +14,12 @@
 
 #include <CGAL/license/AABB_tree.h>
 
-
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Aff_transformation_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Cartesian/Aff_transformation_3.h>
+#include <CGAL/Cartesian_converter.h>
+#include <CGAL/Interval_nt.h>
 #include <CGAL/AABB_tree/internal/Has_nested_type_Shared_data.h>
 #include <boost/mpl/has_xxx.hpp>
 
@@ -61,6 +66,53 @@ struct Primitive_helper<AABBTraits,false>{
   static Reference_point_type get_reference_point(const typename AABBTraits::Primitive& p,const AABBTraits&) {return p.reference_point();}
 };
 
+
+template<class Kernel>
+Bbox_3 compute_transformed_bbox(const CGAL::Aff_transformation_repC3<Kernel>& at, const Bbox_3& bbox, bool has_rotation)
+{
+  typedef Simple_cartesian<Interval_nt<false>> AK;
+  typedef Cartesian_converter<Kernel, AK>    C2F;
+  C2F c2f;
+
+  AK::Aff_transformation_3 a_at = c2f(CGAL::Aff_transformation_3<Kernel>(at));
+  AK::FT xtrm[6] = { c2f((bbox.min)(0)), c2f((bbox.max)(0)),
+                     c2f((bbox.min)(1)), c2f((bbox.max)(1)),
+                     c2f((bbox.min)(2)), c2f((bbox.max)(2)) };
+
+  if(!has_rotation){
+    AK::Point_3 ps[2];
+    ps[0] = a_at( AK::Point_3(xtrm[0], xtrm[2], xtrm[4]) );
+    ps[1] = a_at( AK::Point_3(xtrm[1], xtrm[3], xtrm[5]) );
+
+    return bbox_3(ps, ps+2);
+  }
+
+  AK::Point_3 ps[8];
+  ps[0] = a_at( AK::Point_3(xtrm[0], xtrm[2], xtrm[4]) );
+  ps[1] = a_at( AK::Point_3(xtrm[0], xtrm[2], xtrm[5]) );
+  ps[2] = a_at( AK::Point_3(xtrm[0], xtrm[3], xtrm[4]) );
+  ps[3] = a_at( AK::Point_3(xtrm[0], xtrm[3], xtrm[5]) );
+
+  ps[4] = a_at( AK::Point_3(xtrm[1], xtrm[2], xtrm[4]) );
+  ps[5] = a_at( AK::Point_3(xtrm[1], xtrm[2], xtrm[5]) );
+  ps[6] = a_at( AK::Point_3(xtrm[1], xtrm[3], xtrm[4]) );
+  ps[7] = a_at( AK::Point_3(xtrm[1], xtrm[3], xtrm[5]) );
+
+  return bbox_3(ps, ps+8);
+}
+
+template<class Kernel>
+bool has_rotation(const CGAL::Aff_transformation_3<Kernel>& at){
+  return  (   at.m(0,1) != 0 || at.m(0,2) != 0 || at.m(1,0) != 0
+           || at.m(1,2) != 0 || at.m(2,0) != 0 || at.m(2,1) != 0);
+}
+
+template<class Kernel>
+bool has_rotation(const CGAL::Aff_transformation_repC3<Kernel>& at){
+  return  (   !is_zero(at.cartesian(0,1)) || !is_zero(at.cartesian(0,2)) || !is_zero(at.cartesian(1,0))
+           || !is_zero(at.cartesian(1,2)) || !is_zero(at.cartesian(2,0)) || !is_zero(at.cartesian(2,1)));
+}
+
 } } //namespace CGAL::internal
 
 #endif //CGAL_INTERNAL_AABB_TREE_PRIMITIVE_HELPER
diff --git a/AABB_tree/include/CGAL/AABB_trees/intersection.h b/AABB_tree/include/CGAL/AABB_trees/intersection.h
new file mode 100644
index 000000000000..af2124c54e0e
--- /dev/null
+++ b/AABB_tree/include/CGAL/AABB_trees/intersection.h
@@ -0,0 +1,182 @@
+// Copyright (c) 2026  GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s) : Léo Valque
+
+#ifndef CGAL_AABB_TREES_INTERSECTIONS_H
+#define CGAL_AABB_TREES_INTERSECTIONS_H
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/AABB_tree/internal/AABB_two_trees_traversal.h>
+#include <CGAL/AABB_tree/internal/AABB_two_trees_traversal_traits.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/tbb.h>
+#endif
+
+/// \file AABB_trees/intersection.h
+
+namespace CGAL{
+namespace AABB_trees {
+  /// \ingroup PkgAABBTreeRef
+  ///
+  /// \brief Tests whether two AABB trees contain intersecting primitives.
+  ///
+  /// \cgalNamedParamsBegin
+  ///   \cgalParamNBegin{concurrency_tag}
+  ///     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+  ///     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+  ///     \cgalParamDefault{`CGAL::Sequential_tag`}
+  ///     \cgalParamExtra{`np1` only}
+  ///   \cgalParamNEnd
+  ///   \cgalParamNBegin{transformation}
+  ///     \cgalParamDescription{An affine transformation apply to `tree1` (`tree2`)}
+  ///     \cgalParamType{`CGAL::Aff_transformation_3<Kernel>` where `Kernel` is the kernel associated with `AABBTree1::AABB_traits::Point` (`AABBTree2::AABB_traits::Point`)}
+  ///     \cgalParamDefault{An identity transformation}
+  ///   \cgalParamNEnd
+  /// \cgalNamedParamsEnd
+  ///
+  /// \return `true` if at least one primitive of `tree1` intersects
+  /// a primitive of `tree2`, and `false` otherwise.
+  template< typename AABBTree1,
+            typename AABBTree2,
+            typename NamedParameters1 = parameters::Default_named_parameters,
+            typename NamedParameters2 = parameters::Default_named_parameters>
+  bool do_intersect(const AABBTree1 &tree1,
+                    const AABBTree2 &tree2,
+                    const NamedParameters1& np1 = parameters::default_values(),
+                    const NamedParameters2& np2 = parameters::default_values())
+  {
+    using parameters::get_parameter;
+    using parameters::choose_parameter;
+    using parameters::is_default_parameter;
+    using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                          internal_np::concurrency_tag_t,
+                                          NamedParameters1,
+                                          Sequential_tag
+                                        > ::type;
+
+    if constexpr(is_default_parameter<NamedParameters1, internal_np::transformation_t>::value &&
+                 is_default_parameter<NamedParameters2, internal_np::transformation_t>::value)
+    {
+      CGAL::internal::AABB_tree::Two_trees_do_intersect_traits traversal_traits(tree1.traits(), tree2.traits());
+      CGAL::internal::AABB_tree::two_trees_traversal<Concurrency_tag>(tree1, tree2, traversal_traits);
+      return traversal_traits.is_intersection_found();
+    }
+    else
+    {
+      using Kernel = typename Kernel_traits<typename AABBTree1::AABBTraits::Point>::Kernel::type;
+      using Aff_tr = Aff_transformation_3<Kernel>;
+      const Aff_tr& tr1 = choose_parameter(get_parameter(np1, internal_np::transformation), Aff_tr(Identity_transformation()));
+      const Aff_tr& tr2 = choose_parameter(get_parameter(np2, internal_np::transformation), Aff_tr(Identity_transformation()));
+      CGAL::internal::AABB_tree::Two_trees_do_intersect_traits_with_transformation traversal_traits(tree1.traits(), tree2.traits(), tr1, tr2);
+      // TODO Aff_transformation is not thread safe, no concurrency_tag available
+      CGAL::internal::AABB_tree::two_trees_traversal(tree1, tree2, traversal_traits);
+      return traversal_traits.is_intersection_found();
+    }
+  }
+
+  /// \ingroup PkgAABBTreeRef
+  ///
+  /// \brief Computes all intersecting primitive pairs between two AABB trees.
+  ///
+  /// Traverses both trees and outputs all pairs of primitives that intersect.
+  /// Each output element is a pair `(id1, id2)` where:
+  ///       - `id1` is the ID of a primitive from `tree1`
+  ///       - `id2` is the ID of a primitive from `tree2`
+  ///
+  /// \tparam AABBTree1        Type of the first AABB tree.
+  /// \tparam AABBTree2        Type of the second AABB tree.
+  /// \tparam OutputIterator   Output iterator storing pairs of primitive IDs.
+  /// \tparam NamedParameters1 a sequence of \ref bgl_namedparameters "Named Parameters"
+  /// \tparam NamedParameters2 a sequence of \ref bgl_namedparameters "Named Parameters"
+  ///
+  /// \cgalNamedParamsBegin
+  ///   \cgalParamNBegin{concurrency_tag}
+  ///     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+  ///     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+  ///     \cgalParamDefault{`CGAL::Sequential_tag`}
+  ///     \cgalParamExtra{`np1` only}
+  ///   \cgalParamNEnd
+  ///   \cgalParamNBegin{transformation}
+  ///     \cgalParamDescription{An affine transformation apply to `tree1` (`tree2`)}
+  ///     \cgalParamType{`CGAL::Aff_transformation_3<Kernel>` where `Kernel` is the kernel associated with `AABBTree1::AABB_traits::Point` (`AABBTree2::AABB_traits::Point`)}
+  ///     \cgalParamDefault{An identity transformation}
+  ///   \cgalParamNEnd
+  /// \cgalNamedParamsEnd
+  ///
+  /// \see do_intersect()
+  template< typename AABBTree1,
+            typename AABBTree2,
+            typename OutputIterator,
+            typename NamedParameters1 = parameters::Default_named_parameters,
+            typename NamedParameters2 = parameters::Default_named_parameters>
+  void all_pairs_of_intersecting_primitives(const AABBTree1 &tree1,
+                                            const AABBTree2 &tree2,
+                                            OutputIterator out,
+                                            const NamedParameters1& np1 = parameters::default_values(),
+                                            const NamedParameters2& np2 = parameters::default_values())
+  {
+    using parameters::get_parameter;
+    using parameters::choose_parameter;
+    using parameters::is_default_parameter;
+    using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                          internal_np::concurrency_tag_t,
+                                          NamedParameters1,
+                                          Sequential_tag
+                                        > ::type;
+
+    if constexpr(is_default_parameter<NamedParameters1, internal_np::transformation_t>::value &&
+                 is_default_parameter<NamedParameters2, internal_np::transformation_t>::value)
+    {
+      CGAL::internal::AABB_tree::Two_trees_listing_intersecting_primitives_traits traversal_traits(tree1.traits(), tree2.traits(), out);
+      CGAL::internal::AABB_tree::two_trees_traversal<Concurrency_tag>(tree1, tree2, traversal_traits);
+    }
+    else
+    {
+      using Kernel = typename Kernel_traits<typename AABBTree1::AABB_traits::Point>::Kernel;
+      using Aff_tr = Aff_transformation_3<Kernel>;
+      using Aff_tr_rep = Aff_transformation_repC3<Kernel>;
+      Aff_tr tr1_ref = choose_parameter(get_parameter(np1, internal_np::transformation), Aff_tr(Identity_transformation()));
+      Aff_tr tr2_ref = choose_parameter(get_parameter(np2, internal_np::transformation), Aff_tr(Identity_transformation()));
+      Aff_tr_rep tr1(tr1_ref.m(0,0), tr1_ref.m(0,1), tr1_ref.m(0,2), tr1_ref.m(0,3),
+                     tr1_ref.m(1,0), tr1_ref.m(1,1), tr1_ref.m(1,2), tr1_ref.m(1,3),
+                     tr1_ref.m(2,0), tr1_ref.m(2,1), tr1_ref.m(2,2), tr1_ref.m(2,3));
+      Aff_tr_rep tr2(tr2_ref.m(0,0), tr2_ref.m(0,1), tr2_ref.m(0,2), tr2_ref.m(0,3),
+                     tr2_ref.m(1,0), tr2_ref.m(1,1), tr2_ref.m(1,2), tr2_ref.m(1,3),
+                     tr2_ref.m(2,0), tr2_ref.m(2,1), tr2_ref.m(2,2), tr2_ref.m(2,3));
+      CGAL::internal::AABB_tree::Two_trees_listing_intersecting_primitives_traits_with_transformation traversal_traits(tree1.traits(), tree2.traits(), out, tr1, tr2);
+      // TODO Aff_transformation is not thread safe, no concurrency_tag available
+      CGAL::internal::AABB_tree::two_trees_traversal<Concurrency_tag>(tree1, tree2, traversal_traits);
+    }
+  }
+
+  /// \ingroup PkgAABBTreeRef
+  ///
+  /// \brief Computes all self-intersecting primitive pairs in a single AABB tree.
+  ///
+  /// Intersections of a primitive with itself are not reported, and each intersecting
+  /// pair of distinct primitives is reported only once.
+  template< typename Concurrency_tag = Sequential_tag,
+            typename AABBTree,
+            typename OutputIterator>
+  void all_pairs_of_intersecting_primitives(const AABBTree &tree,
+                                            OutputIterator out)
+  {
+    CGAL::internal::AABB_tree::Listing_self_intersecting_primitives_traits traversal_traits(tree.traits(), out);
+    CGAL::internal::AABB_tree::two_trees_traversal<Concurrency_tag>(tree, tree, traversal_traits);
+  }
+
+}} // end namespace CGAL::AABB_trees
+
+#endif
diff --git a/AABB_tree/test/AABB_tree/CMakeLists.txt b/AABB_tree/test/AABB_tree/CMakeLists.txt
index 1e13ff0ba265..750fc7bf2ddf 100644
--- a/AABB_tree/test/AABB_tree/CMakeLists.txt
+++ b/AABB_tree/test/AABB_tree/CMakeLists.txt
@@ -14,3 +14,11 @@ file(
 foreach(cppfile ${cppfiles})
   create_single_source_cgal_program("${cppfile}")
 endforeach()
+
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(aabb_test_two_trees_intersection PRIVATE CGAL::TBB_support)
+else()
+  message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.")
+endif()
\ No newline at end of file
diff --git a/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp b/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp
new file mode 100644
index 000000000000..8fddeac3f9ad
--- /dev/null
+++ b/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp
@@ -0,0 +1,45 @@
+#include <fstream>
+#include <iostream>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_trees/intersection.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/AABB_tree.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+void test(const std::string fname1, const std::string fname2, std::size_t nb_inter)
+{
+  typedef CGAL::Surface_mesh<K::Point_3> Surface_mesh;
+  typedef CGAL::AABB_face_graph_triangle_primitive<Surface_mesh> Primitive;
+  typedef CGAL::AABB_traits_3<K, Primitive> Traits;
+  typedef CGAL::AABB_tree<Traits> Tree;
+
+  Surface_mesh sm1, sm2;
+  CGAL::IO::read_polygon_mesh(fname1, sm1);
+  CGAL::IO::read_polygon_mesh(fname2, sm2);
+
+  Tree tree1(faces(sm1).begin(), faces(sm1).end(), sm1);
+  Tree tree2(faces(sm2).begin(), faces(sm2).end(), sm2);
+
+  std::vector< std::pair<Primitive::Id, Primitive::Id> > inter;
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter));
+  assert( CGAL::AABB_trees::do_intersect(tree1, tree2) == (nb_inter!=0) );
+  assert( inter.size() == nb_inter );
+}
+
+int main()
+{
+    test(CGAL::data_file_path("meshes/cube.off"), CGAL::data_file_path("meshes/cylinder.off"), 121);
+    test(CGAL::data_file_path("meshes/cube.off"), CGAL::data_file_path("meshes/femur.off"), 0);
+    test(CGAL::data_file_path("meshes/cube.off"), CGAL::data_file_path("meshes/pinion_small.off"), 0);
+    test(CGAL::data_file_path("meshes/cylinder.off"), CGAL::data_file_path("meshes/femur.off"), 0);
+    test(CGAL::data_file_path("meshes/cylinder.off"), CGAL::data_file_path("meshes/pinion_small.off"), 0);
+    test(CGAL::data_file_path("meshes/femur.off"), CGAL::data_file_path("meshes/pinion_small.off"), 905);
+    return EXIT_SUCCESS;
+}
diff --git a/AABB_tree/test/AABB_tree/aabb_test_triangle_soup.cpp b/AABB_tree/test/AABB_tree/aabb_test_triangle_soup.cpp
new file mode 100644
index 000000000000..c74cc809328c
--- /dev/null
+++ b/AABB_tree/test/AABB_tree/aabb_test_triangle_soup.cpp
@@ -0,0 +1,87 @@
+#include <iostream>
+#include <vector>
+#include <array>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/AABB_indexed_triangle_primitive_3.h>
+
+#include <CGAL/AABB_trees/intersection.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 P;
+typedef K::Line_3 L;
+typedef K::Triangle_3 T;
+typedef K::Ray_3 R;
+
+typedef std::vector<std::array<std::size_t, 3> >::const_iterator IndexIterator;
+typedef std::vector<P> PointRange;
+typedef CGAL::AABB_indexed_triangle_primitive_3<K, IndexIterator, PointRange> Primitive;
+typedef CGAL::AABB_traits_3<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+typedef Tree::Point_and_primitive_id Point_and_primitive_id;
+
+int main()
+{
+  P a(0.0, 0.0, 0.0);
+  P b(0.0, 1.0, 0.0);
+  P c(1.0, 0.0, 0.0);
+  P d(1.0, 1.0, 0.0);
+  P e(2.0, 0.0, 0.0);
+  P f(2.0, 1.0, 0.0);
+
+  std::vector<P> points = { a, b, c, d, e, f };
+
+  std::vector<std::array<std::size_t, 3> > triangles;
+  triangles.push_back({ 0, 2, 1 });
+  triangles.push_back({ 1, 2, 3 });
+  triangles.push_back({ 3, 2, 4 });
+  triangles.push_back({ 3, 4, 5 });
+
+  // constructs AABB tree
+  Tree tree(triangles.begin(), triangles.end(), points);
+
+  // point sampling
+  Point_and_primitive_id id;
+  id = tree.closest_point_and_primitive(P(0.5, 0.4, 0));
+  assert(std::distance(triangles.cbegin(), id.second) == 0);
+  id = tree.closest_point_and_primitive(P(0.5, 0.6, 0));
+  assert(std::distance(triangles.cbegin(), id.second) == 1);
+  id = tree.closest_point_and_primitive(P(1.5, 0.4, 0));
+  assert(std::distance(triangles.cbegin(), id.second) == 2);
+  id = tree.closest_point_and_primitive(P(1.5, 0.6, 0));
+  assert(std::distance(triangles.cbegin(), id.second) == 3);
+  id = tree.closest_point_and_primitive(P(3.0, 0.5, 0));
+  assert(std::distance(triangles.cbegin(), id.second) == 3);
+
+  R ray(P(5.5, 0.5, 0), P(1.5, 0.4, 0));
+  auto intersection = tree.first_intersection(ray);
+
+  assert(intersection.has_value());
+  assert(std::distance(triangles.cbegin(), intersection->second) == 3);
+
+  std::vector<P> pts1, pts2;
+  std::vector<std::array<std::size_t, 3> > trs1, trs2;
+  if(!CGAL::IO::read_polygon_soup(CGAL::data_file_path("meshes/knot1.off"), pts1, trs1)){
+    std::cout << "error reading knot1" << std::endl;
+    exit(1);
+  }
+  if(!CGAL::IO::read_polygon_soup(CGAL::data_file_path("meshes/lion.off"), pts2, trs2)){
+    std::cout << "error reading lion" << std::endl;
+    exit(1);
+  }
+
+  Tree tree1(trs1.begin(), trs1.end(), pts1);
+  Tree tree2(trs2.begin(), trs2.end(), pts2);
+  tree1.build();
+  tree2.build();
+
+  using Iterator = std::vector<std::array<std::size_t, 3> >::const_iterator;
+  std::vector< std::pair<Iterator, Iterator> > inter;
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter));
+  assert(inter.size() == 1191);
+
+  return EXIT_SUCCESS;
+}
diff --git a/AABB_tree/test/AABB_tree/aabb_test_two_trees_intersection.cpp b/AABB_tree/test/AABB_tree/aabb_test_two_trees_intersection.cpp
new file mode 100644
index 000000000000..273251d3373d
--- /dev/null
+++ b/AABB_tree/test/AABB_tree/aabb_test_two_trees_intersection.cpp
@@ -0,0 +1,113 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Simple_cartesian.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+
+#include <CGAL/AABB_trees/intersection.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+
+using Epick = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Epeck = CGAL::Exact_predicates_exact_constructions_kernel;
+using Cartesian = CGAL::Simple_cartesian<double>;
+
+template<typename K>
+void test()
+{
+  using P = typename K::Point_3;
+  using V = typename K::Vector_3;
+  using M =  CGAL::Surface_mesh<P>;
+  using fd = typename boost::graph_traits<M>::face_descriptor;
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<M>;
+  using Traits = CGAL::AABB_traits_3<K, Primitive>;
+  using Tree = CGAL::AABB_tree<Traits>;
+  using Aff_tr = CGAL::Aff_transformation_3<K>;
+
+  M m1, m2;
+  if(!CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/knot1.off"), m1)){
+    std::cout << "error reading knot1" << std::endl;
+    exit(1);
+  }
+  if(!CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/lion.off"), m2)){
+    std::cout << "error reading lion" << std::endl;
+    exit(1);
+  }
+
+  Tree tree1(faces(m1).first, faces(m1).second, m1);
+  Tree tree2(faces(m2).first, faces(m2).second, m2);
+  tree1.build();
+  tree2.build();
+
+  std::vector< std::pair<fd, fd> > inter;
+  assert(CGAL::AABB_trees::do_intersect(tree1, tree2));
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter));
+  assert(inter.size() == 1191);
+  inter.clear();
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), CGAL::parameters::transformation(Aff_tr(CGAL::Translation(), V(1,0,0))));
+  assert(inter.size() == 0);
+  inter.clear();
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), CGAL::parameters::transformation(Aff_tr(0, 1, 0, 1, 0, 0, 0, 0, 1, 1)));
+  assert(inter.size() == 1289);
+}
+
+#ifdef CGAL_LINKED_WITH_TBB
+template<typename K>
+void test_parallel()
+{
+  using P = typename K::Point_3;
+  using V = typename K::Vector_3;
+  using M =  CGAL::Surface_mesh<P>;
+  using fd = typename boost::graph_traits<M>::face_descriptor;
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<M>;
+  using Traits = CGAL::AABB_traits_3<K, Primitive>;
+  using Tree = CGAL::AABB_tree<Traits>;
+  using Aff_tr = CGAL::Aff_transformation_3<K>;
+
+  M m1, m2;
+  if(!CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/knot1.off"), m1)){
+    std::cout << "error reading knot1" << std::endl;
+    exit(1);
+  }
+  if(!CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/lion.off"), m2)){
+    std::cout << "error reading lion" << std::endl;
+    exit(1);
+  }
+
+  Tree tree1(faces(m1).first, faces(m1).second, m1);
+  Tree tree2(faces(m2).first, faces(m2).second, m2);
+  tree1.template build<CGAL::Parallel_tag>();
+  tree2.template build<CGAL::Parallel_tag>();
+
+  tbb::concurrent_vector< std::pair<fd, fd> > inter;
+  assert(CGAL::AABB_trees::do_intersect(tree1, tree2, CGAL::parameters::concurrency_tag(CGAL::Parallel_tag())));
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()));
+  std::cout << inter.size() << std::endl;
+  assert(inter.size() == 1191);
+  inter.clear();
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()).transformation(Aff_tr(CGAL::Translation(), V(0.5,0,0))));
+  std::cout << inter.size() << std::endl;
+  assert(inter.size() == 0);
+  inter.clear();
+
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives(tree1, tree2, std::back_inserter(inter), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()).transformation(Aff_tr(0, 1, 0, 1, 0, 0, 0, 0, 1, 1)));
+  std::cout << inter.size() << std::endl;
+  assert(inter.size() == 1289);
+}
+#endif
+
+int main(){
+  test<Epick>();
+  test<Epeck>();
+  test<Cartesian>();
+
+  return 0;
+}
diff --git a/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_3.h b/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_3.h
index d696e1917a30..bb6005c11d2a 100644
--- a/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_3.h
+++ b/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_3.h
@@ -40,7 +40,7 @@ template < class R_ >
 class Aff_transformationC3
   : public Handle_for_virtual<Aff_transformation_rep_baseC3<R_> >
 {
-  friend class PlaneC3<R_>; // FIXME: why ?
+  // friend class PlaneC3<R_>; // FIXME: why ?
 
   typedef typename R_::FT                   FT;
   typedef Aff_transformation_rep_baseC3<R_> Aff_t_base;
@@ -117,6 +117,12 @@ class Aff_transformationC3
                                             m31, m32, m33, m34));
   }
 
+  // General form: with translation
+  Aff_transformationC3(const Aff_transformation_repC3<R> &b)
+  {
+    initialize_with(b);
+  }
+
   Point_3
   transform(const Point_3 &p) const
   { return this->Ptr()->transform(p); }
diff --git a/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_rep_3.h b/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_rep_3.h
index 2e441e194d3b..9c272d90a83a 100644
--- a/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_rep_3.h
+++ b/Cartesian_kernel/include/CGAL/Cartesian/Aff_transformation_rep_3.h
@@ -135,7 +135,7 @@ class Aff_transformation_repC3
 
   virtual Plane_3 transform(const Plane_3& p) const
   {
-      if (is_even())
+    if (is_even())
       return Plane_3(transform(p.point()),
                  transpose().inverse().transform(p.orthogonal_direction()));
     else
@@ -157,8 +157,8 @@ class Aff_transformation_repC3
   virtual bool is_even() const
   {
     return sign_of_determinant(t11, t12, t13,
-                                  t21, t22, t23,
-                                  t31, t32, t33) == POSITIVE;
+                               t21, t22, t23,
+                               t31, t32, t33) == POSITIVE;
   }
 
 
diff --git a/Cartesian_kernel/include/CGAL/Cartesian/function_objects.h b/Cartesian_kernel/include/CGAL/Cartesian/function_objects.h
index 56390484622d..cc30818c03e2 100644
--- a/Cartesian_kernel/include/CGAL/Cartesian/function_objects.h
+++ b/Cartesian_kernel/include/CGAL/Cartesian/function_objects.h
@@ -4122,6 +4122,17 @@ namespace CartesianKernelFunctors {
     typedef typename K::Sphere_3       Sphere_3;
 
   public:
+    std::pair<Orientation, Orientation>
+    operator()( const Point_3& p, const Point_3& q,
+                const Point_3& r, const Point_3& s, const Point_3& t) const
+    {
+      return orientationsC3(p.x(), p.y(), p.z(),
+                            q.x(), q.y(), q.z(),
+                            r.x(), r.y(), r.z(),
+                            s.x(), s.y(), s.z(),
+                            t.x(), t.y(), t.z());
+    }
+
     Orientation
     operator()( const Point_3& p, const Point_3& q,
                 const Point_3& r, const Point_3& s) const
diff --git a/Cartesian_kernel/include/CGAL/predicates/kernel_ftC3.h b/Cartesian_kernel/include/CGAL/predicates/kernel_ftC3.h
index 2d6fe0e53be1..8ac90b350fa1 100644
--- a/Cartesian_kernel/include/CGAL/predicates/kernel_ftC3.h
+++ b/Cartesian_kernel/include/CGAL/predicates/kernel_ftC3.h
@@ -112,6 +112,21 @@ collinearC3(const FT &px, const FT &py, const FT &pz,
                                  sign_of_determinant(dpy, dqy, dpz, dqz) == ZERO );
 }
 
+
+template < class FT >
+CGAL_KERNEL_MEDIUM_INLINE
+typename Same_uncertainty_nt<std::pair<Orientation,Orientation>, FT>::type
+orientationsC3(const FT &px, const FT &py, const FT &pz,
+               const FT &qx, const FT &qy, const FT &qz,
+               const FT &rx, const FT &ry, const FT &rz,
+               const FT &sx, const FT &sy, const FT &sz,
+               const FT &tx, const FT &ty, const FT &tz)
+{
+  return sign_of_determinants<FT>(qx-px,rx-px,sx-px,tx-px,
+                                  qy-py,ry-py,sy-py,ty-py,
+                                  qz-pz,rz-pz,sz-pz,tz-pz);
+}
+
 template < class FT >
 CGAL_KERNEL_MEDIUM_INLINE
 typename Same_uncertainty_nt<Orientation, FT>::type
@@ -121,8 +136,8 @@ orientationC3(const FT &px, const FT &py, const FT &pz,
               const FT &sx, const FT &sy, const FT &sz)
 {
   return sign_of_determinant<FT>(qx-px,rx-px,sx-px,
-                                    qy-py,ry-py,sy-py,
-                                    qz-pz,rz-pz,sz-pz);
+                                 qy-py,ry-py,sy-py,
+                                 qz-pz,rz-pz,sz-pz);
 }
 
 template < class FT >
diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h
index ce73ed3e83c5..c34a45fea492 100644
--- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h
+++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h
@@ -41,6 +41,138 @@ class Orientation_3
 public:
   using Base::operator();
 
+  std::pair<Orientation,Orientation>
+  operator()(const Point_3 &p, const Point_3 &q,
+             const Point_3 &r, const Point_3 &s, const Point_3 &t) const
+  {
+      CGAL_BRANCH_PROFILER_3("semi-static failures/attempts/calls to   : Orientation_3", tmp);
+
+      double px, py, pz, qx, qy, qz, rx, ry, rz, sx, sy, sz, tx, ty, tz;
+
+      if (fit_in_double(p.x(), px) && fit_in_double(p.y(), py) &&
+          fit_in_double(p.z(), pz) &&
+          fit_in_double(q.x(), qx) && fit_in_double(q.y(), qy) &&
+          fit_in_double(q.z(), qz) &&
+          fit_in_double(r.x(), rx) && fit_in_double(r.y(), ry) &&
+          fit_in_double(r.z(), rz) &&
+          fit_in_double(s.x(), sx) && fit_in_double(s.y(), sy) &&
+          fit_in_double(s.z(), sz) &&
+          fit_in_double(t.x(), tx) && fit_in_double(t.y(), ty) &&
+          fit_in_double(t.z(), tz))
+      {
+          CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
+
+          double pqx = qx - px;
+          double pqy = qy - py;
+          double pqz = qz - pz;
+          double prx = rx - px;
+          double pry = ry - py;
+          double prz = rz - pz;
+          double psx = sx - px;
+          double psy = sy - py;
+          double psz = sz - pz;
+          double ptx = tx - px;
+          double pty = ty - py;
+          double ptz = tz - pz;
+
+          // CGAL::abs uses fabs on platforms where it is faster than (a<0)?-a:a
+          // Then semi-static filter.
+
+          double maxx = CGAL::abs(pqx);
+          double maxy = CGAL::abs(pqy);
+          double maxz = CGAL::abs(pqz);
+
+          double aprx = CGAL::abs(prx);
+          double apsx = CGAL::abs(psx);
+          double aptx = CGAL::abs(ptx);
+
+          double apry = CGAL::abs(pry);
+          double apsy = CGAL::abs(psy);
+          double apty = CGAL::abs(pty);
+
+          double aprz = CGAL::abs(prz);
+          double apsz = CGAL::abs(psz);
+          double aptz = CGAL::abs(ptz);
+#ifdef CGAL_USE_SSE2_MAX
+          CGAL::Max<double> mmax;
+
+          maxx = mmax(maxx, aprx, apsx); // todo add aptx etc
+          maxy = mmax(maxy, apry, apsy);
+          maxz = mmax(maxz, aprz, apsz);
+#else
+
+
+          if (maxx < aprx) maxx = aprx;
+          double maxxt = maxx;
+          if (maxx < apsx) maxx = apsx;
+          if (maxxt < aptx) maxxt = aptx;
+          if (maxy < apry) maxy = apry;
+          double maxyt = maxy;
+          if (maxy < apsy) maxy = apsy;
+          if (maxyt < apty) maxyt= apty;
+          if (maxz < aprz) maxz = aprz;
+          double maxzt = maxz;
+          if (maxz < apsz) maxz = apsz;
+          if (maxzt < aptz) maxzt = aptz;
+#endif
+          std::pair<double,double> det = CGAL::determinants(pqx, prx, psx, ptx,
+                                                            pqy, pry, psy, pty,
+                                                            pqz, prz, psz, ptz);
+
+          double epss = 5.1107127829973299e-15 * maxx * maxy * maxz;
+          double epst = 5.1107127829973299e-15 * maxxt * maxyt * maxzt;
+
+#ifdef CGAL_USE_SSE2_MAX
+#if 0
+          CGAL::Min<double> mmin;
+          double tmp = mmin(maxx, maxy, maxz);
+          maxz = mmax(maxx, maxy, maxz);
+          maxx = tmp;
+#else
+          sse2minmax(maxx,maxy,maxz);
+          // maxy can contain ANY element
+#endif
+#else
+          // Sort maxx < maxy < maxz.
+          if (maxx > maxz)
+              std::swap(maxx, maxz);
+          if (maxy > maxz)
+              std::swap(maxy, maxz);
+          else if (maxy < maxx)
+              std::swap(maxx, maxy);
+#endif
+          std::pair<Orientation,Orientation> res = std::make_pair(ZERO,ZERO);
+          bool first = false;
+          bool second = false;
+          // Protect against underflow in the computation of eps.
+          if (maxx < 1e-97) /* cbrt(min_double/eps) */ {
+            if (maxx == 0)
+              first = true;
+          }
+          if (maxxt < 1e-97) /* cbrt(min_double/eps) */ {
+            if (maxxt == 0)
+              second = true;
+          }
+          // Protect against overflow in the computation of det.
+          if ((!first) && (maxz < 1e102)) /* cbrt(max_double [hadamard]/4) */ {
+             if (det.first > epss)  { res.first = POSITIVE; first = true;}
+             else if (det.first < -epss) { res.first = NEGATIVE; first = true;}
+          }
+
+          if((! second) && (maxzt < 1e102)){
+            if (det.second > epst) { res.second = POSITIVE; second = true;}
+               else if (det.second < -epst) { res.second = NEGATIVE; second = true;}
+          }
+          if(first && second){
+            return res;
+          }
+          }
+          CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
+    return Base::operator()(p, q, r, s, t);
+  }
+
+
+
   Orientation
   operator()(const Point_3 &p, const Point_3 &q,
              const Point_3 &r, const Point_3 &s) const
diff --git a/Homogeneous_kernel/include/CGAL/Homogeneous/function_objects.h b/Homogeneous_kernel/include/CGAL/Homogeneous/function_objects.h
index 6c14ddd06b51..a353deb8fbb8 100644
--- a/Homogeneous_kernel/include/CGAL/Homogeneous/function_objects.h
+++ b/Homogeneous_kernel/include/CGAL/Homogeneous/function_objects.h
@@ -4347,6 +4347,14 @@ namespace HomogeneousKernelFunctors {
     typedef typename K::Sphere_3       Sphere_3;
 
   public:
+
+    std::pair<Orientation, Orientation>
+    operator()( const Point_3& p, const Point_3& q,
+                const Point_3& r, const Point_3& s, const Point_3& t) const
+    {
+      return std::make_pair( (*this)(p, q, r, s), (*this)(p, q, r, t) );
+    }
+
     Orientation
     operator()( const Point_3& p, const Point_3& q,
                 const Point_3& r, const Point_3& s) const
diff --git a/Intersections_3/benchmark/Intersections_3/CMakeLists.txt b/Intersections_3/benchmark/Intersections_3/CMakeLists.txt
new file mode 100644
index 000000000000..9b37b7281dda
--- /dev/null
+++ b/Intersections_3/benchmark/Intersections_3/CMakeLists.txt
@@ -0,0 +1,16 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.12...3.31)
+project(Intersections_3_Benchmark)
+
+find_package(CGAL REQUIRED)
+
+# create a target per cppfile
+file(
+  GLOB cppfiles
+  RELATIVE ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+foreach(cppfile ${cppfiles})
+  create_single_source_cgal_program("${cppfile}")
+endforeach()
diff --git a/Intersections_3/benchmark/Intersections_3/segment_triangle_3.cpp b/Intersections_3/benchmark/Intersections_3/segment_triangle_3.cpp
new file mode 100644
index 000000000000..cd355dcbebaa
--- /dev/null
+++ b/Intersections_3/benchmark/Intersections_3/segment_triangle_3.cpp
@@ -0,0 +1,50 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Timer.h>
+#include <iostream>
+
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = Kernel::Point_3;
+using Segment_3 = Kernel::Segment_3;
+using Triangle_3 = Kernel::Triangle_3;
+
+int main()
+{
+  Point_3 p(0, 0, 0);
+  Point_3 q(10, 0, 0);
+  Point_3 r(0, 10, 0);
+
+  Point_3 s(1, 1, -1);
+  Point_3 t(1, 1, 1);
+  Point_3 u(1, 1, -2);
+
+  Segment_3 st(s, t);
+  Segment_3 tu(t, u);
+  Triangle_3 tri(p, q, r);
+
+  int does = 0;
+  int does_not = 0;
+
+   CGAL::Timer timer;
+   timer.start();
+
+
+  for (int i = 0; i < 1000000000; ++i) {
+    int j = i/3;
+    if(i == 3*j){
+      if(CGAL::do_intersect(st, tri)){
+        ++does;
+      }
+    }else{
+      if(CGAL::do_intersect(tu, tri)){
+        ++does_not;
+      }
+    }
+  }
+
+   timer.stop();
+   std::cout << "Does intersect: " << does << std::endl;
+   std::cout << "Does not intersect: " << does_not << std::endl;
+   std::cout << "Time: " << timer.time() << " seconds." << std::endl;
+
+  return 0;
+}
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Segment_3_Triangle_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Segment_3_Triangle_3_do_intersect.h
index e595c0d01a8e..9a8198bcc1e3 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Segment_3_Triangle_3_do_intersect.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Segment_3_Triangle_3_do_intersect.h
@@ -194,8 +194,19 @@ do_intersect(const typename K::Triangle_3& t,
   const Point_3& p = point_on(s,0);
   const Point_3& q = point_on(s,1);
 
-  const Orientation abcp = orientation(a,b,c,p);
-  const Orientation abcq = orientation(a,b,c,q);
+#if 1
+  std::pair<Orientation,Orientation> abcp_abcq = orientation(a,b,c,p,q);
+
+  const Orientation abcp = abcp_abcq.first;
+  const Orientation abcq = abcp_abcq.second;
+  CGAL_assertion_code(const Orientation abcpbis = orientation(a,b,c,p);)
+  CGAL_assertion_code(const Orientation abcqbis = orientation(a,b,c,q);)
+  CGAL_assertion(abcpbis == abcp);
+  CGAL_assertion(abcqbis == abcq);
+#else
+const Orientation abcp = orientation(a,b,c,p);
+const Orientation abcq = orientation(a,b,c,q);
+#endif
 
   switch ( abcp ) {
   case POSITIVE:
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h b/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
index fb0de9bc83eb..804b7a9109f2 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
@@ -87,6 +87,26 @@ double ymax() const;
 */
 double zmax() const;
 
+/*!
+
+*/
+double x_span() const;
+
+/*!
+
+*/
+double y_span() const;
+
+/*!
+
+*/
+double z_span() const;
+
+/*!
+
+*/
+double squared_diagonal_length() const;
+
 /*!
 Returns `xmin()` if `i==0` or `ymin()` if `i==1`
 or `zmin()` if `i==2`.
diff --git a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
index 27e8a539d8e4..c4f5c22e2cd1 100644
--- a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
+++ b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
@@ -9442,6 +9442,15 @@ class Orientation_3 {
                          const Kernel::Point_3&r,
                          const Kernel::Point_3&s);
 
+  /*!
+    performs the same test as the previous operator simultaneously for `(p,q,r,s)` and `(p,q,r,t)`.
+  */
+  std::pair<Orientation,Orientation> operator()(const Kernel::Point_3&p,
+                                                const Kernel::Point_3&q,
+                                                const Kernel::Point_3&r,
+                                                const Kernel::Point_3&s,
+                                                const Kernel::Point_3&t);
+
   /*!
     returns \ref CGAL::POSITIVE if `u`, `v` and `w` are positively oriented,
     returns \ref CGAL::NEGATIVE if `u`, `v` and `w` are negatively oriented,
diff --git a/Kernel_23/include/CGAL/Bbox_3.h b/Kernel_23/include/CGAL/Bbox_3.h
index dc750e4d0e0e..e1c9f22051ae 100644
--- a/Kernel_23/include/CGAL/Bbox_3.h
+++ b/Kernel_23/include/CGAL/Bbox_3.h
@@ -66,6 +66,7 @@ class Bbox_3
   double x_span() const;
   double y_span() const;
   double z_span() const;
+  double squared_diagonal_length() const;
 
   inline double min BOOST_PREVENT_MACRO_SUBSTITUTION (int i) const;
   inline double max BOOST_PREVENT_MACRO_SUBSTITUTION (int i) const;
@@ -122,6 +123,10 @@ inline double Bbox_3::z_span() const {
   return zmax() - zmin();
 }
 
+inline double Bbox_3::squared_diagonal_length() const {
+  return x_span()*x_span() + y_span()*y_span() + z_span()*z_span();
+}
+
 inline
 bool
 Bbox_3::operator==(const Bbox_3 &b) const
diff --git a/Kernel_23/include/CGAL/Kernel/Same_uncertainty.h b/Kernel_23/include/CGAL/Kernel/Same_uncertainty.h
index db32f878f970..f3398c92492b 100644
--- a/Kernel_23/include/CGAL/Kernel/Same_uncertainty.h
+++ b/Kernel_23/include/CGAL/Kernel/Same_uncertainty.h
@@ -36,6 +36,12 @@ struct Same_uncertainty < T1, Uncertain<T2> >
   typedef Uncertain<T1> type;
 };
 
+template < typename T1, typename T2 >
+struct Same_uncertainty < std::pair<T1,T1>, Uncertain<T2> >
+{
+  typedef std::pair<Uncertain<T1>,Uncertain<T1> > type;
+};
+
 // Short cut to extract uncertainty from a number type directly.
 template < typename T1, typename NT >
 struct Same_uncertainty_nt
diff --git a/Kernel_23/include/CGAL/determinant.h b/Kernel_23/include/CGAL/determinant.h
index f470720d5609..b39d69337f94 100644
--- a/Kernel_23/include/CGAL/determinant.h
+++ b/Kernel_23/include/CGAL/determinant.h
@@ -19,6 +19,7 @@
 #define CGAL_DETERMINANT_H
 
 #include <CGAL/kernel_config.h>
+#include <utility>
 
 namespace CGAL {
 
@@ -41,7 +42,7 @@ determinant(
  const RT& a00,  const RT& a01,  const RT& a02,
  const RT& a10,  const RT& a11,  const RT& a12,
  const RT& a20,  const RT& a21,  const RT& a22)
-{
+ {
 // First compute the det2x2
   const RT m01 = a00*a11 - a10*a01;
   const RT m02 = a00*a21 - a20*a01;
@@ -51,6 +52,25 @@ determinant(
   return m012;
 }
 
+template <class RT>
+CGAL_KERNEL_MEDIUM_INLINE
+std::pair<RT,RT>
+determinants(
+ const RT& a00,  const RT& a01,  const RT& a02, const RT& b02,
+ const RT& a10,  const RT& a11,  const RT& a12, const RT& b12,
+ const RT& a20,  const RT& a21,  const RT& a22, const RT& b22)
+{
+// First compute the det2x2
+  const RT m01 = a00*a11 - a10*a01;
+  const RT m02 = a00*a21 - a20*a01;
+  const RT m12 = a10*a21 - a20*a11;
+// Now compute the minors of rank 3
+  const RT m012a = m01*a22 - m02*a12 + m12*a02;
+  const RT m012b = m01*b22 - m02*b12 + m12*b02;
+  return std::make_pair(m012a, m012b);
+}
+
+
 template <class RT>
 CGAL_KERNEL_LARGE_INLINE
 RT
diff --git a/Kernel_23/include/CGAL/predicates/sign_of_determinant.h b/Kernel_23/include/CGAL/predicates/sign_of_determinant.h
index 5543772cfb3a..0ecc4db095db 100644
--- a/Kernel_23/include/CGAL/predicates/sign_of_determinant.h
+++ b/Kernel_23/include/CGAL/predicates/sign_of_determinant.h
@@ -46,6 +46,21 @@ sign_of_determinant( const RT& a00,  const RT& a01,  const RT& a02,
                                          a20, a21, a22));
 }
 
+template <class RT>
+inline
+std::pair<typename Sgn<RT>::result_type,typename Sgn<RT>::result_type>
+sign_of_determinants(const RT& a00, const RT& a01, const RT& a02, const RT& b02,
+                     const RT& a10, const RT& a11, const RT& a12, const RT& b12,
+                     const RT& a20, const RT& a21, const RT& a22, const RT& b22)
+{
+  std::pair<RT,RT> det = determinants(a00, a01, a02, b02,
+                                      a10, a11, a12, b12,
+                                      a20, a21, a22, b22);
+
+ return std::make_pair(CGAL_NTS sign(det.first), CGAL_NTS sign(det.second));
+}
+
+
 template <class RT>
 inline
 typename Sgn<RT>::result_type
diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt
index ede30624b931..f3ebdd2b83ab 100644
--- a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/CMakeLists.txt
@@ -42,3 +42,12 @@ else()
 endif()
 
 create_single_source_cgal_program("fast.cpp")
+
+create_single_source_cgal_program("meshes_intersections_comparison.cpp")
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(meshes_intersections_comparison PRIVATE CGAL::TBB_support)
+else()
+  message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.")
+endif()
diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/include/AABB_meshes_intersections.h b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/include/AABB_meshes_intersections.h
new file mode 100644
index 000000000000..79762ccf7b17
--- /dev/null
+++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/include/AABB_meshes_intersections.h
@@ -0,0 +1,737 @@
+// Copyright (c) 2026  GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s) : Léo Valque
+
+#ifndef CGAL_AABB_MESHES_INTERSECTIONS_H
+#define CGAL_AABB_MESHES_INTERSECTIONS_H
+
+#include <memory>
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+
+#include <CGAL/AABB_trees/intersection.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/Box_intersection_d/Box_with_info_d.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/tbb.h>
+#include <CGAL/mutex.h>
+#endif
+
+namespace CGAL {
+
+namespace Polygon_mesh_processing{
+
+namespace experimental{
+
+template< typename TriangleMesh1,
+          typename TriangleMesh2,
+          typename OutputIterator,
+          typename NamedParameters1 = parameters::Default_named_parameters,
+          typename NamedParameters2 = parameters::Default_named_parameters>
+void mixed_meshes_intersections(const TriangleMesh1 &tm1,
+                                        const TriangleMesh2 &tm2,
+                                        OutputIterator out,
+                                        const NamedParameters1& np1 = parameters::default_values(),
+                                        const NamedParameters2& np2 = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
+  using halfedge_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::halfedge_descriptor;
+
+  using face_descriptor_2 = typename boost::graph_traits<TriangleMesh2>::face_descriptor;
+  using halfedge_descriptor_2 = typename boost::graph_traits<TriangleMesh2>::halfedge_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters1,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh1, NamedParameters1>::type;
+  GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+
+  using Face_bbox_tag = typename CGAL::dynamic_face_property_t<Bbox_3>;
+  using Primitive_1 = AABB_face_graph_triangle_primitive<TriangleMesh1>;
+  using Bbox_pmap_1 = typename boost::property_map<TriangleMesh1, Face_bbox_tag>::const_type;
+  using Traits_1 = AABB_traits_3<GT, Primitive_1, Bbox_pmap_1>;
+  using Tree_1 = AABB_tree<Traits_1>;
+  using Node_1 = AABB_node<Traits_1>;
+  using ConstPrimitiveIterator_1 = typename std::vector< Primitive_1 >::const_iterator;
+  using Box_1 = Box_intersection_d::Box_with_info_d<double, 3, face_descriptor_1>;
+
+  using Primitive_2 = AABB_face_graph_triangle_primitive<TriangleMesh2>;
+  using Bbox_pmap_2 = typename boost::property_map<TriangleMesh2, Face_bbox_tag>::const_type;
+  using Traits_2 = AABB_traits_3<GT, Primitive_2, Bbox_pmap_2>;
+  using Tree_2 = AABB_tree<Traits_2>;
+  using Node_2 = AABB_node<Traits_2>;
+  using ConstPrimitiveIterator_2 = typename std::vector< Primitive_2 >::const_iterator;
+  using Box_2 = Box_intersection_d::Box_with_info_d<double, 3, face_descriptor_2>;
+
+  using InternOutputIterator= std::back_insert_iterator<std::vector<std::pair<Node_1*, Node_2*>>>;
+
+  const std::size_t cutoff = 50000;
+
+  auto vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm1));
+  auto vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm2));
+
+  auto triangle = [&](auto fd, const auto &vpm, const auto &tm){
+    auto hd = halfedge(fd,tm);
+    auto a = get(vpm, source(hd,tm));
+    auto b = get(vpm, target(hd,tm));
+    auto c = get(vpm, target(next(hd,tm),tm));
+    return typename GT::Triangle_3(a, b, c);
+  };
+
+  auto bbox = [](auto fd, const auto &vpm, const auto &tm){
+    auto hd = halfedge(fd,tm);
+    Bbox_3 res = get(vpm, source(hd,tm)).bbox();
+    res += get(vpm, target(hd,tm)).bbox();
+    res += get(vpm, target(next(hd,tm),tm)).bbox();
+    return res;
+  };
+
+  Bbox_pmap_1 bb1 = get(Face_bbox_tag(), tm1);
+  Bbox_pmap_2 bb2 = get(Face_bbox_tag(), tm2);
+#ifdef CGAL_LINKED_WITH_TBB
+  if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+  {
+    oneapi::tbb::parallel_for(
+      oneapi::tbb::blocked_range<size_t>(0, faces(tm1).size()),
+        [&](const oneapi::tbb::blocked_range<size_t>& r) {
+          for (size_t i = r.begin(); i < r.end(); ++i) {
+            face_descriptor_1 fd = *(faces(tm1).begin() + i);
+            put(bb1, fd, bbox(fd, vpm1, tm1));
+          }
+        }
+    );
+
+    oneapi::tbb::parallel_for(
+      oneapi::tbb::blocked_range<size_t>(0, faces(tm2).size()),
+        [&](const oneapi::tbb::blocked_range<size_t>& r) {
+          for (size_t i = r.begin(); i < r.end(); ++i) {
+            face_descriptor_2 fd = *(faces(tm2).begin() + i);
+            put(bb2, fd, bbox(fd, vpm2, tm2));
+          }
+        }
+    );
+  }
+  else
+#endif
+  {
+    for(face_descriptor_1 fd : faces(tm1))
+      put(bb1, fd, bbox(fd, vpm1, tm1));
+    for(face_descriptor_2 fd : faces(tm2))
+      put(bb2, fd, bbox(fd, vpm2, tm2));
+  }
+
+  Traits_1 traits1(bb1);
+  Tree_1 tree1(traits1);
+  tree1.insert(faces(tm1).first, faces(tm1).second, tm1);
+
+  Traits_2 traits2(bb2);
+  Tree_2 tree2(traits2);
+  tree2.insert(faces(tm2).first, faces(tm2).second, tm2);
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+  {
+    CGAL_MUTEX m;
+    oneapi::tbb::task_group tg;
+    tg.run([&]{ tree1.template partial_build<Concurrency_tag>(cutoff); });
+    tree2.template partial_build<Concurrency_tag>(cutoff);
+    tg.wait();
+
+    tbb::concurrent_vector<std::pair<const Node_1*, const Node_2*>> inter;
+    CGAL::internal::AABB_tree::experimental::Two_trees_intersecting_nodes_traits traversal_traits(tree1.traits(), tree2.traits(), std::back_inserter(inter));
+    CGAL::internal::AABB_tree::experimental::two_trees_partial_traversal<Concurrency_tag>(tree1, tree2, cutoff, traversal_traits);
+    oneapi::tbb::parallel_for(
+      oneapi::tbb::blocked_range<size_t>(0, inter.size()),
+      [&](const oneapi::tbb::blocked_range<size_t>& r) {
+        for (size_t i = r.begin(); i != r.end(); ++i) {
+          const auto [begin_1, end_1] = tree1.partial_node_to_primitives_iterator(*inter[i].first);
+          const auto [begin_2, end_2] = tree2.partial_node_to_primitives_iterator(*inter[i].second);
+          std::vector< Box_1 > boxes_1;
+          std::vector< Box_2 > boxes_2;
+          std::vector< Box_1* > boxes_ptr_1;
+          std::vector< Box_2* > boxes_ptr_2;
+          boxes_1.reserve(std::distance(begin_1, end_1));
+          boxes_2.reserve(std::distance(begin_2, end_2));
+          boxes_ptr_1.reserve(boxes_1.size());
+          boxes_ptr_2.reserve(boxes_2.size());
+
+          for(auto it=begin_1; it!=end_1; ++it)
+            boxes_1.emplace_back(get(bb1, it->id()), it->id());
+          for(auto it=begin_2; it!=end_2; ++it)
+            boxes_2.emplace_back(get(bb2, it->id()), it->id());
+          for(auto &b: boxes_1)
+            boxes_ptr_1.push_back(std::addressof(b));
+          for(auto &b: boxes_2)
+            boxes_ptr_2.push_back(std::addressof(b));
+
+          std::vector< std::pair<face_descriptor_1, face_descriptor_2>> patch_out;
+          const std::ptrdiff_t cutoff = 2000;
+          box_intersection_d(boxes_ptr_1.begin(), boxes_ptr_1.end(), boxes_ptr_2.begin(), boxes_ptr_2.end(),
+              [&](const Box_1 *b1, const Box_2 *b2)
+              {
+                if(CGAL::do_intersect(triangle(b1->info(), vpm1, tm1), triangle(b2->info(), vpm2, tm2)))
+                  patch_out.emplace_back(b1->info(), b2->info());
+              }, cutoff);
+              CGAL_SCOPED_LOCK(m);
+          for(auto p: patch_out)
+            *out ++ = p;
+        }
+      }
+    );
+  }
+  else
+#endif
+  {
+    tree1.template partial_build<Concurrency_tag>(cutoff);
+    tree2.template partial_build<Concurrency_tag>(cutoff);
+
+    std::vector<std::pair<const Node_1*, const Node_2*>> inter;
+    CGAL::internal::AABB_tree::experimental::Two_trees_intersecting_nodes_traits traversal_traits(tree1.traits(), tree2.traits(), std::back_inserter(inter));
+    CGAL::internal::AABB_tree::experimental::two_trees_partial_traversal(tree1, tree2, cutoff, traversal_traits);
+
+    for(const auto& [n_1, n_2]: inter){
+      const auto [begin_1, end_1] = tree1.partial_node_to_primitives_iterator(*n_1);
+      const auto [begin_2, end_2] = tree2.partial_node_to_primitives_iterator(*n_2);
+      std::vector< Box_1 > boxes_1;
+      std::vector< Box_2 > boxes_2;
+      std::vector< Box_1* > boxes_ptr_1;
+      std::vector< Box_2* > boxes_ptr_2;
+      boxes_1.reserve(std::distance(begin_1, end_1));
+      boxes_2.reserve(std::distance(begin_2, end_2));
+      boxes_ptr_1.reserve(boxes_1.size());
+      boxes_ptr_2.reserve(boxes_2.size());
+
+      for(auto it=begin_1; it!=end_1; ++it)
+        boxes_1.emplace_back(get(bb1, it->id()), it->id());
+      for(auto it=begin_2; it!=end_2; ++it)
+        boxes_2.emplace_back(get(bb2, it->id()), it->id());
+      for(auto &b: boxes_1)
+        boxes_ptr_1.push_back(std::addressof(b));
+      for(auto &b: boxes_2)
+        boxes_ptr_2.push_back(std::addressof(b));
+
+      std::vector< std::pair<face_descriptor_1, face_descriptor_2>> patch_out;
+      const std::ptrdiff_t cutoff = 2000;
+      box_intersection_d(boxes_ptr_1.begin(), boxes_ptr_1.end(), boxes_ptr_2.begin(), boxes_ptr_2.end(),
+          [&](const Box_1 *b1, const Box_2 *b2)
+          {
+            if(CGAL::do_intersect(triangle(b1->info(), vpm1, tm1), triangle(b2->info(), vpm2, tm2)))
+              patch_out.emplace_back(b1->info(), b2->info());
+          }, cutoff);
+      for(auto p: patch_out)
+        *out ++ = p;
+    }
+  }
+}
+
+template< typename TriangleMesh1,
+          typename TriangleMesh2,
+          typename OutputIterator,
+          typename NamedParameters1 = parameters::Default_named_parameters,
+          typename NamedParameters2 = parameters::Default_named_parameters>
+void box_meshes_intersections(const TriangleMesh1 &tm1,
+                                        const TriangleMesh2 &tm2,
+                                        OutputIterator out,
+                                        const NamedParameters1& np1 = parameters::default_values(),
+                                        const NamedParameters2& np2 = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
+  using halfedge_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::halfedge_descriptor;
+
+  using face_descriptor_2 = typename boost::graph_traits<TriangleMesh2>::face_descriptor;
+  using halfedge_descriptor_2 = typename boost::graph_traits<TriangleMesh2>::halfedge_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters1,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh1, NamedParameters1>::type;
+  GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+
+  using Box_1 = Box_intersection_d::Box_with_info_d<double, 3, face_descriptor_1>;
+  using Box_2 = Box_intersection_d::Box_with_info_d<double, 3, face_descriptor_2>;
+
+  auto vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm1));
+  auto vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm2));
+
+  auto triangle = [&](auto fd, const auto &vpm, const auto &tm){
+    auto hd = halfedge(fd,tm);
+    auto a = get(vpm, source(hd,tm));
+    auto b = get(vpm, target(hd,tm));
+    auto c = get(vpm, target(next(hd,tm),tm));
+    return typename GT::Triangle_3(a, b, c);
+  };
+
+  auto bbox = [](auto fd, const auto &vpm, const auto &tm){
+    auto hd = halfedge(fd,tm);
+    Bbox_3 res = get(vpm, source(hd,tm)).bbox();
+    res += get(vpm, target(hd,tm)).bbox();
+    res += get(vpm, target(next(hd,tm),tm)).bbox();
+    return res;
+  };
+
+  std::vector< Box_1 > boxes_1;
+  std::vector< Box_2 > boxes_2;
+  boxes_1.reserve(faces(tm1).size());
+  boxes_2.reserve(faces(tm2).size());
+  for(auto fd: faces(tm1))
+    boxes_1.emplace_back(bbox(fd, vpm1, tm1), fd);
+  for(auto fd: faces(tm2))
+    boxes_2.emplace_back(bbox(fd, vpm2, tm2), fd);
+
+  std::vector< Box_1* > boxes_ptr_1;
+  std::vector< Box_2* > boxes_ptr_2;
+  boxes_ptr_1.reserve(faces(tm1).size());
+  boxes_ptr_2.reserve(faces(tm2).size());
+  for(auto &b: boxes_1)
+    boxes_ptr_1.emplace_back(std::addressof(b));
+  for(auto &b: boxes_2)
+    boxes_ptr_2.emplace_back(std::addressof(b));
+
+  CGAL_MUTEX m;
+  auto callback=[&](const Box_1 *b1, const Box_2 *b2){
+    if(CGAL::do_intersect(triangle(b1->info(), vpm1, tm1), triangle(b2->info(), vpm2, tm2))){
+      CGAL_SCOPED_LOCK(m);
+      *out ++ = std::make_pair(b1->info(), b2->info());
+    }
+  };
+
+  const std::ptrdiff_t cutoff = 2000;
+  box_intersection_d<Concurrency_tag>(boxes_ptr_1.begin(), boxes_ptr_1.end(), boxes_ptr_2.begin(), boxes_ptr_2.end(), callback, cutoff);
+}
+
+template< typename TriangleMesh,
+          typename OutputIterator,
+          typename NamedParameters = parameters::Default_named_parameters>
+void AABB_two_trees_self_intersections(const TriangleMesh &tm, OutputIterator out, const NamedParameters& np = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh, NamedParameters>::type;
+  GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+  using Primitive = AABB_face_graph_triangle_primitive<TriangleMesh>;
+  using Face_bbox_tag = typename CGAL::dynamic_face_property_t<Bbox_3>  ;
+  using Bbox_pmap = typename boost::property_map<TriangleMesh, Face_bbox_tag>::const_type;
+  using Traits = AABB_traits_3<GT, Primitive, Bbox_pmap>;
+  using Tree = AABB_tree<Traits>;
+
+  auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                              get_const_property_map(vertex_point, tm));
+
+  auto bbox = [&](face_descriptor fd){
+    halfedge_descriptor hd = halfedge(fd,tm);
+    Bbox_3 res = get(vpm, source(hd,tm)).bbox();
+    res += get(vpm, target(hd,tm)).bbox();
+    res += get(vpm, target(next(hd,tm),tm)).bbox();
+    return res;
+  };
+
+  Bbox_pmap bb = get(Face_bbox_tag(), tm);
+  for(face_descriptor fd : faces(tm))
+    put(bb, fd, bbox(fd));
+
+  Traits traits(bb);
+  Tree tree(traits);
+  tree.insert(faces(tm).first, faces(tm).second, tm);
+  tree.template build<Concurrency_tag>();
+
+#ifdef CGAL_LINKED_WITH_TBB
+
+#endif
+  using InternOutputIterator= std::back_insert_iterator<std::vector<std::pair<face_descriptor, face_descriptor>>>;
+  std::vector<std::pair<face_descriptor, face_descriptor>> inter;
+  CGAL::AABB_trees::all_pairs_of_intersecting_primitives<Concurrency_tag>(tree, std::back_inserter(inter));
+  for(const auto& [f_1, f_2]: inter)
+    if(f_1 < f_2)
+      if(Polygon_mesh_processing::internal::do_faces_intersect<GT>(f_1, f_2, tm, tm.points(), gt.construct_segment_3_object(), gt.construct_triangle_3_object(), gt.do_intersect_3_object()))
+        *out ++ = std::make_pair(f_1, f_2);
+}
+
+template <class RPM>
+struct Split_primitives
+{
+  Split_primitives(RPM rpm)
+    : rpm(rpm)
+  {}
+
+  template<typename PrimitiveIterator>
+  void operator()(PrimitiveIterator first,
+                  PrimitiveIterator beyond,
+                  const CGAL::Bbox_3& bbox) const
+  {
+    auto longest_axis=[](const CGAL::Bbox_3& bbox){
+      const double dx = bbox.xmax() - bbox.xmin();
+      const double dy = bbox.ymax() - bbox.ymin();
+      const double dz = bbox.zmax() - bbox.zmin();
+      return (dx>=dy) ? ((dx>=dz) ? 0 : 2) : ((dy>=dz) ? 1 : 2);
+    };
+
+    PrimitiveIterator middle = first + (beyond - first)/2;
+    typedef typename std::iterator_traits<PrimitiveIterator>::value_type Primitive;
+    const int crd=longest_axis(bbox);
+    const RPM& l_rpm=rpm;
+    std::nth_element(first, middle, beyond,
+                    [l_rpm, crd](const Primitive& p1, const Primitive& p2){ return get(l_rpm, p1.id())[crd] < get(l_rpm, p2.id())[crd];});
+  }
+  RPM rpm;
+};
+
+template <class BBM>
+struct Compute_bbox {
+  Compute_bbox(const BBM& bbm)
+    : bbm(bbm)
+  {}
+
+  template<typename ConstPrimitiveIterator>
+  CGAL::Bbox_3 operator()(ConstPrimitiveIterator first,
+                          ConstPrimitiveIterator beyond) const
+  {
+    CGAL::Bbox_3 bbox = get(bbm, first->id());
+    for(++first; first != beyond; ++first)
+    {
+      bbox += get(bbm, first->id());
+    }
+    return bbox;
+  }
+  BBM bbm;
+};
+
+template< typename TriangleMesh1,
+          typename TriangleMesh2,
+          typename OutputIterator,
+          typename NamedParameters1 = parameters::Default_named_parameters,
+          typename NamedParameters2 = parameters::Default_named_parameters>
+void AABB_two_trees_meshes_intersections(const TriangleMesh1 &tm1,
+                                        const TriangleMesh2 &tm2,
+                                        OutputIterator out,
+                                        const NamedParameters1& np1 = parameters::default_values(),
+                                        const NamedParameters2& np2 = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
+  using face_descriptor_2 = typename boost::graph_traits<TriangleMesh2>::face_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters1,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh1, NamedParameters1>::type;
+  // GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+
+  using Primitive_1 = AABB_face_graph_triangle_primitive<TriangleMesh1>;
+  using Traits_1 = AABB_traits_3<GT, Primitive_1>;
+  using Tree_1 = AABB_tree<Traits_1>;
+
+  using Primitive_2 = AABB_face_graph_triangle_primitive<TriangleMesh2>;
+  using Traits_2 = AABB_traits_3<GT, Primitive_2>;
+  using Tree_2 = AABB_tree<Traits_2>;
+
+  auto vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm1));
+  auto vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm2));
+
+    // Custom build functor using pointers
+    // TODO THIS IS SURFACE MESH SPECIFIC
+
+    typedef Pointer_property_map<Bbox_3>::type BBM;
+    typedef Pointer_property_map<Epick::Point_3>::type RPM; // EPIC on purpose here
+
+    std::vector<Bbox_3> v_bb1;
+    std::vector<Bbox_3> v_bb2;
+    std::vector<Epick::Point_3> v_rp1;
+    std::vector<Epick::Point_3> v_rp2;
+
+    std::size_t nbf1 = faces(tm1).size();
+    std::size_t nbf2 = faces(tm2).size();
+    v_bb1.resize(nbf1);
+    v_bb2.resize(nbf2);
+    v_rp1.resize(nbf1);
+    v_rp2.resize(nbf2);
+    BBM bbmap1 = make_property_map(v_bb1);
+    BBM bbmap2 = make_property_map(v_bb2);
+    RPM rpm1 = make_property_map(v_rp1);
+    RPM rpm2 = make_property_map(v_rp2);
+
+    CGAL::Cartesian_converter<GT, Epick> to_input;
+  #ifdef CGAL_LINKED_WITH_TBB
+    if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+    {
+      tbb::parallel_for(std::size_t(0), faces(tm1).size(), [&](std::size_t i){
+        face_descriptor_1 f(i);
+        v_bb1[f]=face_bbox(f, tm1);
+        v_rp1[f]=to_input(get(vpm1, (target(halfedge(f, tm1), tm1))));
+      });
+      tbb::parallel_for(std::size_t(0), faces(tm2).size(), [&](std::size_t i){
+        face_descriptor_2 f(i);
+        v_bb2[f]=face_bbox(f, tm2);
+        v_rp2[f]=to_input(get(vpm2, (target(halfedge(f, tm2), tm2))));
+      });
+    }
+    else
+  #endif
+    {
+      for(face_descriptor_1 f : faces(tm1)){
+        v_bb1[f]=face_bbox(f, tm1);
+        v_rp1[f]=to_input(get(vpm1, (target(halfedge(f, tm1), tm1))));
+      }
+      for(face_descriptor_2 f : faces(tm2)){
+        v_bb2[f]=face_bbox(f, tm2);
+        v_rp2[f]=to_input(get(vpm2, (target(halfedge(f, tm2), tm2))));
+      }
+    }
+
+    Compute_bbox<BBM> compute_bbox1(bbmap1);
+    Compute_bbox<BBM> compute_bbox2(bbmap2);
+    Split_primitives<RPM> split_primitives1(rpm1);
+    Split_primitives<RPM> split_primitives2(rpm2);
+
+    Tree_1 tree1;
+    Tree_2 tree2;
+    tree1.insert(faces(tm1).first, faces(tm1).second, tm1);
+    tree2.insert(faces(tm2).first, faces(tm2).second, tm2);
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+  {
+    oneapi::tbb::task_group tg;
+    tg.run([&]{ tree1.template custom_build<Concurrency_tag>(compute_bbox1, split_primitives1); });
+    tree2.template custom_build<Concurrency_tag>(compute_bbox2, split_primitives2);
+    tg.wait();
+
+    CGAL::AABB_trees::all_pairs_of_intersecting_primitives<Concurrency_tag>(tree1, tree2, out);
+  }
+  else
+#endif
+  {
+    tree1.template custom_build(compute_bbox1, split_primitives1);
+    tree2.template custom_build(compute_bbox2, split_primitives2);
+    CGAL::AABB_trees::all_pairs_of_intersecting_primitives<Concurrency_tag>(tree1, tree2, out);
+  }
+}
+
+} // end of namespace experimental
+
+template< typename TriangleMesh,
+          typename OutputIterator,
+          typename NamedParameters = parameters::Default_named_parameters>
+void AABB_self_intersections(const TriangleMesh &tm, OutputIterator out, const NamedParameters& np = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh, NamedParameters>::type;
+  GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+  using Primitive = AABB_face_graph_triangle_primitive<TriangleMesh>;
+  using Face_bbox_tag = typename CGAL::dynamic_face_property_t<Bbox_3>  ;
+  using Bbox_pmap = typename boost::property_map<TriangleMesh, Face_bbox_tag>::const_type;
+  using Traits = AABB_traits_3<GT, Primitive, Bbox_pmap>;
+  using Tree = AABB_tree<Traits>;
+
+  auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                              get_const_property_map(vertex_point, tm));
+
+  auto bbox = [&](face_descriptor fd){
+    halfedge_descriptor hd = halfedge(fd,tm);
+    Bbox_3 res = get(vpm, source(hd,tm)).bbox();
+    res += get(vpm, target(hd,tm)).bbox();
+    res += get(vpm, target(next(hd,tm),tm)).bbox();
+    return res;
+  };
+
+  Bbox_pmap bb = get(Face_bbox_tag(), tm);
+  for(face_descriptor fd : faces(tm))
+    put(bb, fd, bbox(fd));
+
+  Traits traits(bb);
+  Tree tree(traits);
+  tree.insert(faces(tm).first, faces(tm).second, tm);
+  tree.template build<Concurrency_tag>();
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+  {
+    CGAL_MUTEX mutex;
+    tbb::concurrent_vector<std::pair<face_descriptor, face_descriptor>> inter;
+    std::vector<face_descriptor> face_vec(faces(tm).begin(), faces(tm).end());
+    oneapi::tbb::parallel_for(
+      oneapi::tbb::blocked_range<std::size_t>(0, face_vec.size()),
+      [&](const oneapi::tbb::blocked_range<std::size_t>& r) {
+        for (std::size_t i = r.begin(); i != r.end(); ++i) {
+          face_descriptor f_1 = face_vec[i];
+          std::vector<face_descriptor> inter;
+          tree.all_intersected_primitives(get(bb, f_1), std::back_inserter(inter));
+          for(face_descriptor f_2: inter)
+            if(f_1 < f_2)
+              if(Polygon_mesh_processing::internal::do_faces_intersect<GT>(f_1, f_2, tm, tm.points(), gt.construct_segment_3_object(), gt.construct_triangle_3_object(), gt.do_intersect_3_object())){
+                CGAL_SCOPED_LOCK(mutex);
+                *out ++ = std::make_pair(f_1, f_2);
+              }
+        }
+      }
+    );
+  }
+  else
+#endif
+  {
+    for (face_descriptor f_1: faces(tm)){
+      std::vector<face_descriptor> inter;
+      tree.all_intersected_primitives(get(bb, f_1), std::back_inserter(inter));
+      for(face_descriptor f_2: inter)
+        if(f_1 < f_2)
+          if(Polygon_mesh_processing::internal::do_faces_intersect<GT>(f_1, f_2, tm, tm.points(), gt.construct_segment_3_object(), gt.construct_triangle_3_object(), gt.do_intersect_3_object()))
+            *out ++ = std::make_pair(f_1, f_2);
+    }
+  }
+}
+
+template< typename TriangleMesh1,
+          typename TriangleMesh2,
+          typename OutputIterator,
+          typename NamedParameters1 = parameters::Default_named_parameters,
+          typename NamedParameters2 = parameters::Default_named_parameters>
+void AABB_meshes_intersections(const TriangleMesh1 &tm1,
+                               const TriangleMesh2 &tm2,
+                               OutputIterator out,
+                               const NamedParameters1& np1 = parameters::default_values(),
+                               const NamedParameters2& np2 = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using face_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
+  using halfedge_descriptor_1 = typename boost::graph_traits<TriangleMesh1>::halfedge_descriptor;
+
+  using face_descriptor_2 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
+  using halfedge_descriptor_2 = typename boost::graph_traits<TriangleMesh1>::halfedge_descriptor;
+
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def <
+                                                internal_np::concurrency_tag_t,
+                                                NamedParameters1,
+                                                Sequential_tag
+                                              > ::type;
+  using GT = typename GetGeomTraits<TriangleMesh1, NamedParameters1>::type;
+  GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+
+  using Primitive = AABB_face_graph_triangle_primitive<TriangleMesh2>;
+  using Face_bbox_tag = typename CGAL::dynamic_face_property_t<Bbox_3>  ;
+  using Bbox_pmap = typename boost::property_map<TriangleMesh2, Face_bbox_tag>::const_type;
+  using Traits = AABB_traits_3<GT, Primitive, Bbox_pmap>;
+  using Tree = AABB_tree<Traits>;
+
+  auto vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm1));
+  auto vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tm2));
+
+  auto triangle = [&](face_descriptor_1 fd){
+    halfedge_descriptor_1 hd = halfedge(fd,tm1);
+    auto a = get(vpm1, source(hd,tm1));
+    auto b = get(vpm1, target(hd,tm1));
+    auto c = get(vpm1, target(next(hd,tm1),tm1));
+    return typename GT::Triangle_3(a, b, c);
+  };
+
+  auto bbox = [&](face_descriptor_2 fd){
+    halfedge_descriptor_2 hd = halfedge(fd,tm2);
+    Bbox_3 res = get(vpm2, source(hd,tm2)).bbox();
+    res += get(vpm2, target(hd,tm2)).bbox();
+    res += get(vpm2, target(next(hd,tm2),tm2)).bbox();
+    return res;
+  };
+
+  Bbox_pmap bb = get(Face_bbox_tag(), tm2);
+  for(face_descriptor_2 fd : faces(tm2))
+    put(bb, fd, bbox(fd));
+
+  Traits traits(bb);
+  Tree tree(traits);
+  tree.insert(faces(tm2).first, faces(tm2).second, tm2);
+  tree.template build<Concurrency_tag>();
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if constexpr(std::is_same_v<Concurrency_tag, Parallel_tag>)
+  {
+    CGAL_MUTEX mutex;
+    tbb::concurrent_vector<std::pair<face_descriptor_1, face_descriptor_2>> inter;
+    std::vector<face_descriptor_1> face_vec(faces(tm1).begin(), faces(tm1).end());
+    oneapi::tbb::parallel_for(
+      oneapi::tbb::blocked_range<size_t>(0, face_vec.size()),
+      [&](const oneapi::tbb::blocked_range<size_t>& r) {
+        for (size_t i = r.begin(); i != r.end(); ++i) {
+          face_descriptor_1 f_1 = face_vec[i];
+          std::vector<face_descriptor_2> inter;
+          tree.all_intersected_primitives(triangle(f_1), std::back_inserter(inter));
+          CGAL_SCOPED_LOCK(mutex);
+          for(face_descriptor_2 f_2: inter)
+            *out ++ = std::make_pair(f_1, f_2);
+        }
+      }
+    );
+  }
+  else
+#endif
+  {
+    for (face_descriptor_1 f_1: faces(tm1)){
+      std::vector<face_descriptor_2> inter;
+      tree.all_intersected_primitives(triangle(f_1), std::back_inserter(inter));
+      for(face_descriptor_2 f_2: inter)
+        *out ++ = std::make_pair(f_1, f_2);
+    }
+  }
+}
+
+}
+} // end namespace CGAL
+
+#endif
diff --git a/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/meshes_intersections_comparison.cpp b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/meshes_intersections_comparison.cpp
new file mode 100644
index 000000000000..98c6e42e0750
--- /dev/null
+++ b/Polygon_mesh_processing/benchmark/Polygon_mesh_processing/meshes_intersections_comparison.cpp
@@ -0,0 +1,148 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/corefinement.h>
+
+#include <CGAL/AABB_tree/internal/AABB_two_trees_traversal.h>
+
+#include <tbb/tbb.h>
+
+#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
+
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/boost/graph/generators.h>
+
+#include <CGAL/Polygon_mesh_processing/transform.h>
+#include <CGAL/Aff_transformation_3.h>
+
+#include "include/AABB_meshes_intersections.h"
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+
+#include <iostream>
+#include <fstream>
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
+
+typedef boost::graph_traits<Surface_mesh>::face_descriptor vertex_descriptor;
+typedef boost::graph_traits<Surface_mesh>::face_descriptor face_descriptor;
+
+void two_meshes_intersection(std::string fname1, std::string fname2){
+
+  Surface_mesh tm1;
+  Surface_mesh tm2;
+  CGAL::IO::read_polygon_mesh(fname1, tm1);
+  CGAL::IO::read_polygon_mesh(fname2, tm2);
+  PMP::triangulate_faces(tm1);
+  PMP::triangulate_faces(tm2);
+
+  CGAL::Timer t;
+  CGAL::Real_timer rt;
+  t.start();
+  rt.start();
+
+#if CGAL_LINKED_WITH_TBB
+  tbb::concurrent_vector<std::pair<face_descriptor, face_descriptor>> out;
+#else
+  std::vector<std::pair<face_descriptor, face_descriptor>> out;
+#endif
+
+  out.clear();
+  PMP::experimental::AABB_two_trees_meshes_intersections(tm1, tm2, std::back_inserter(out), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Two tree AABB intersecton time: " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::experimental::AABB_two_trees_meshes_intersections(tm1, tm2, std::back_inserter(out));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Two tree AABB intersecton time (Sequential): " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::experimental::mixed_meshes_intersections(tm1, tm2, std::back_inserter(out), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Mixed intersecton time: " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::experimental::mixed_meshes_intersections(tm1, tm2, std::back_inserter(out));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Mixed intersecton time (Sequential): " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::AABB_meshes_intersections(tm1, tm2, std::back_inserter(out), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "AABB intersecton time: " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+  std::sort(out.begin(), out.end(), [](const auto &a, const auto &b){ return a.first<b.first || (a.first == b.first && a.second<b.second); });
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::AABB_meshes_intersections(tm1, tm2, std::back_inserter(out));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "AABB intersecton time (Sequential): " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::experimental::box_meshes_intersections(tm1, tm2, std::back_inserter(out), CGAL::parameters::concurrency_tag(CGAL::Parallel_tag()));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Box intersecton time: " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+  std::sort(out.begin(), out.end(), [](const auto &a, const auto &b){ return a.first<b.first || (a.first == b.first && a.second<b.second); });
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+
+  out.clear();
+  PMP::experimental::box_meshes_intersections(tm1, tm2, std::back_inserter(out));
+  std::cout << "number intersections: " << out.size() << std::endl;
+  std::cout << "Box intersecton time (Sequential): " << rt.time() << "sec (" << t.time() << "s all cpu)." << std::endl;
+
+  t.stop(); rt.stop();
+  t.reset(); rt.reset();
+  t.start(); rt.start();
+}
+
+int main(int argc, char* argv[])
+{
+  two_meshes_intersection((argc>1)?argv[1]:CGAL::data_file_path("meshes/tetrahedron.off"),
+                          (argc>2)?argv[2]:CGAL::data_file_path("meshes/beam.off"));
+  // two_meshes_intersection((argc>1)?argv[1]:CGAL::data_file_path("meshes/beam.off"),
+  //                         (argc>2)?argv[2]:"beam_transformed.off");
+  return EXIT_SUCCESS;
+}
diff --git a/Spatial_searching/include/CGAL/Kd_tree.h b/Spatial_searching/include/CGAL/Kd_tree.h
index d23f225713ee..5e8390f10831 100644
--- a/Spatial_searching/include/CGAL/Kd_tree.h
+++ b/Spatial_searching/include/CGAL/Kd_tree.h
@@ -300,11 +300,6 @@ class Kd_tree {
     return pts.empty();
   }
 
-  void build()
-  {
-    build<Sequential_tag>();
-  }
-
   /*
     Note about parallel `build()`. Several different strategies have
     been tried, among which:
@@ -325,7 +320,7 @@ class Kd_tree {
     * the parallel computations are launched using
       `tbb::parallel_invoke`
   */
-  template <typename ConcurrencyTag>
+  template <typename ConcurrencyTag = Sequential_tag>
   void
   build()
   {