diff --git a/AABB_tree/include/CGAL/AABB_tree/internal/triangle_datum_covering.h b/AABB_tree/include/CGAL/AABB_tree/internal/triangle_datum_covering.h
new file mode 100644
index 00000000000..6f7df8ae57c
--- /dev/null
+++ b/AABB_tree/include/CGAL/AABB_tree/internal/triangle_datum_covering.h
@@ -0,0 +1,406 @@
+// Copyright (c) 2021 INRIA Sophia-Antipolis (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) : Mael Rouxel-Labbé
+
+#ifndef CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H
+#define CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+#include <CGAL/AABB_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/array.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Container_helper.h>
+#include <CGAL/property_map.h>
+
+#include <queue>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+namespace CGAL {
+namespace AABB_trees {
+namespace internal {
+
+// std::vector to property map
+// Not using Pointer_property_map because the underlying range is empty initially and can change
+template <typename T>
+struct Vector_property_map
+{
+  using Range = std::vector<T>;
+
+  using key_type = std::size_t;
+  using value_type = T;
+  using reference = value_type&;
+  using category = boost::read_write_property_map_tag;
+
+  Vector_property_map() : m_range_ptr(std::make_shared<Range>()) { }
+
+  inline friend void put(const Vector_property_map& map, const key_type& k, const value_type& v)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+
+    if(k >= map.m_range_ptr->size())
+      map.m_range_ptr->resize(k+1);
+
+    map.m_range_ptr->operator[](k) = v;
+  }
+
+  inline friend reference get(const Vector_property_map& map, const key_type& k)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+    return map.m_range_ptr->operator[](k);
+  }
+
+  Range& range() { return *m_range_ptr; }
+  const Range& range() const { return *m_range_ptr; }
+
+private:
+  std::shared_ptr<Range> m_range_ptr;
+};
+
+// Same as the standard traversal traits, but for multiple primitives per datum,
+// such that the final operation on the datum is only performed once.
+template <typename AABBTraits, typename BaseTraversalTraits>
+struct Covered_traversal_traits
+  : BaseTraversalTraits
+{
+  using Base = BaseTraversalTraits;
+  using Primitive = typename AABBTraits::Primitive;
+
+  std::unordered_set<std::size_t> visited_data;
+
+public:
+  template <typename ... Args>
+  Covered_traversal_traits(Args&&... args) : Base(std::forward<Args>(args)...) { }
+
+  template <typename Query>
+  void intersection(const Query& query, const Primitive& primitive)
+  {
+    // check a datum only once
+    auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
+    if(!is_insert_successful.second)
+      return;
+
+    return Base::intersection(query, primitive);
+  }
+};
+
+template <typename AABBTraits>
+struct Covered_tree_traversal_traits
+{
+  using Base_projection_traits = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
+  using Projection_traits = Covered_traversal_traits<AABBTraits, Base_projection_traits>;
+
+  template <typename Query>
+  using Do_intersect_traits_base = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
+  template <typename Query>
+  using Do_intersect_traits = Covered_traversal_traits<AABBTraits, Do_intersect_traits_base<Query> >;
+
+  template <typename Query>
+  using First_intersection_traits_base = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
+  template <typename Query>
+  using First_intersection_traits = Covered_traversal_traits<AABBTraits, First_intersection_traits_base<Query> >;
+
+  template <typename Query>
+  using First_primitive_traits_base = CGAL::internal::AABB_tree::First_primitive_traits<AABBTraits, Query>;
+  template <typename Query>
+  using First_primitive_traits = Covered_traversal_traits<AABBTraits, First_primitive_traits_base<Query> >;
+
+  template <typename Query, typename CountingIterator>
+  using Listing_primitive_traits_base = CGAL::internal::AABB_tree::Listing_primitive_traits<AABBTraits, Query, CountingIterator>;
+  template <typename Query, typename CountingIterator>
+  using Listing_primitive_traits = Covered_traversal_traits<AABBTraits, Listing_primitive_traits_base<Query, CountingIterator> >;
+
+  template <typename Query, typename CountingIterator>
+  using Listing_intersection_traits_base = CGAL::internal::AABB_tree::Listing_intersection_traits<AABBTraits, Query, CountingIterator>;
+  template <typename Query, typename CountingIterator>
+  using Listing_intersection_traits = Covered_traversal_traits<AABBTraits, Listing_intersection_traits_base<Query, CountingIterator> >;
+};
+
+// Dissociated from the class `AABB_covered_triangle_tree` for clarity
+template <typename Kernel, typename Point>
+struct AABB_covered_triangle_tree_traits
+{
+  using Triangle_3 = typename Kernel::Triangle_3;
+
+  // Below is a lot of trouble to cover a single datum with multiple primitives using smaller bboxes
+  using ID = std::pair<std::size_t /*primitive ID*/, std::size_t /*input face ID*/>;
+  using IDPM = CGAL::First_of_pair_property_map<ID>;
+
+  // Primitive ID --> box vector pos --> Bounding Box
+  using BPMB = internal::Vector_property_map<CGAL::Bbox_3>;
+  using BPM = CGAL::Property_map_binder<IDPM, BPMB>;
+
+  // Primitive ID --> point vector pos --> Reference Point
+  using RPPMB = internal::Vector_property_map<Point>;
+  using RPPM = CGAL::Property_map_binder<IDPM, RPPMB>;
+
+  // Primitive ID --> Datum pos vector pos --> Datum pos --> Datum
+  // The vector of data has size nf, but the vector of datum pos has size tree.size()
+  using DPPMB = internal::Vector_property_map<std::size_t>; // pos --> Datum pos
+  using DPPM = CGAL::Property_map_binder<IDPM, DPPMB>; // PID --> Datum pos
+  using DPMB = internal::Vector_property_map<Triangle_3>; // Datum pos --> Datum
+  using DPM = CGAL::Property_map_binder<DPPM, DPMB>; // PID --> Datum
+
+  using Primitive = CGAL::AABB_primitive<ID, DPM, RPPM,
+                                         CGAL::Tag_true /*external pmaps*/,
+                                         CGAL::Tag_false /*no caching*/>;
+
+  using AABB_geom_traits = Kernel;
+  using AABB_traits = CGAL::AABB_traits<AABB_geom_traits, Primitive, BPM>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+};
+
+template <typename Kernel,
+          typename Point = typename Kernel::Point_3>
+struct AABB_covered_triangle_tree
+  : public AABB_covered_triangle_tree_traits<Kernel, Point>::AABB_tree
+{
+  using FT = typename Kernel::FT;
+  using Triangle_3 = typename Kernel::Triangle_3;
+
+  using ACTTT = AABB_covered_triangle_tree_traits<Kernel, Point>;
+
+  using BPM = typename ACTTT::BPM;
+  using RPPM = typename ACTTT::RPPM;
+  using DPPMB = typename ACTTT::DPPMB;
+  using DPPM = typename ACTTT::DPPM;
+  using DPMB = typename ACTTT::DPMB;
+  using DPM = typename ACTTT::DPM;
+
+  using ID = typename ACTTT::ID;
+  using Primitive = typename ACTTT::Primitive;
+  using AABB_traits = typename ACTTT::AABB_traits;
+  using AABB_tree = typename ACTTT::AABB_tree;
+  using Base = AABB_tree;
+
+protected:
+  double m_sq_length;
+
+  DPPMB m_dppmb; // std::size_t (id) --> datum pos
+
+  BPM m_bpm; // std::size_t (id) --> bounding box
+  RPPM m_rppm; // std::size_t (id) --> reference point
+  DPMB m_dpmb; // std::size_t (datum pos) --> triangle datum
+
+  DPM m_dpm; // std::size_t (id) --> triangle (datum)
+
+  std::size_t fid = 0;
+
+public:
+  AABB_covered_triangle_tree(const double max_length,
+                             const AABB_traits& traits = AABB_traits())
+    : Base(traits),
+      m_sq_length(square(max_length)),
+      m_dppmb(), m_bpm(), m_rppm(), m_dpmb(),
+      m_dpm(DPPM(m_dppmb/*first binder's value_map*/)/*second binder's key map*/, m_dpmb)
+  {
+    initialize_tree_property_maps();
+  }
+
+private:
+  void initialize_tree_property_maps() const
+  {
+    // Can't be set in the default constructed traits that are passed to the base
+    // since m_bpm is a member of the derived class.
+    //
+    // 'const_cast' because CGAL::AABB_tree only gives a const& to its traits...
+    const_cast<AABB_traits&>(this->traits()).bbm = m_bpm;
+    const_cast<AABB_traits&>(this->traits()).set_shared_data(m_dpm, m_rppm);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  CGAL::Bbox_3 compute_bbox(const P& ap0, const P& ap1, const P& ap2)
+  {
+    double xmin = (CGAL::min)(ap0.x().inf(), (CGAL::min)(ap1.x().inf(), ap2.x().inf()));
+    double ymin = (CGAL::min)(ap0.y().inf(), (CGAL::min)(ap1.y().inf(), ap2.y().inf()));
+    double zmin = (CGAL::min)(ap0.z().inf(), (CGAL::min)(ap1.z().inf(), ap2.z().inf()));
+
+    double xmax = (CGAL::max)(ap0.x().sup(), (CGAL::max)(ap1.x().sup(), ap2.x().sup()));
+    double ymax = (CGAL::max)(ap0.y().sup(), (CGAL::max)(ap1.y().sup(), ap2.y().sup()));
+    double zmax = (CGAL::max)(ap0.z().sup(), (CGAL::max)(ap1.z().sup(), ap2.z().sup()));
+
+    return CGAL::Bbox_3(xmin, ymin, zmin, xmax, ymax, zmax);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  const P& compute_reference_point(const P&, const P&, const P& ap2)
+  {
+    return ap2; // ap2 is the midpoint when splitting
+  }
+
+public:
+  void reserve(std::size_t nf)
+  {
+    CGAL::internal::reserve(m_dpmb.range(), m_dpmb.range().size() + nf);
+
+    // Due to splitting, these might need more than 'nf'
+    CGAL::internal::reserve(m_dppmb.range(), m_dppmb.range().size() + nf);
+    CGAL::internal::reserve(m_rppm.value_map.range(), m_rppm.value_map.range().size() + nf);
+    CGAL::internal::reserve(m_bpm.value_map.range(), m_bpm.value_map.range().size() + nf);
+  }
+
+  void split_and_insert(const Triangle_3& tr)
+  {
+    // Convert to intervals to ensure that the bounding box fully covers the subdividing triangle
+    using IT = CGAL::Interval_nt<true>;
+    using NT = typename IT::value_type;
+    using AK = CGAL::Simple_cartesian<IT>;
+    using K2AK = CGAL::Cartesian_converter<Kernel, AK>;
+    using AFT = AK::FT;
+    using APoint_3 = AK::Point_3;
+
+    using APL = std::pair<APoint_3, NT>; // point and upper bound of the length of the opposite edge
+    using AT = std::array<APL, 3>;
+
+    const std::size_t data_size = m_dpmb.range().size();
+    put(m_dpmb, data_size, tr);
+
+    auto vertex = Kernel().construct_vertex_3_object();
+    const Point& p0 = vertex(tr, 0);
+    const Point& p1 = vertex(tr, 1);
+    const Point& p2 = vertex(tr, 2);
+
+    if(m_sq_length == FT(0)) // no splits
+    {
+      const std::size_t pid = this->size();
+      ID id = std::make_pair(pid, fid++);
+
+      put(m_dppmb, pid, data_size);
+      put(m_rppm, id, p1); // the ref point that `One_point_from_face_descriptor_map` would give
+      put(m_bpm, id, Kernel().construct_bbox_3_object()(tr));
+
+//      std::cout << "Primitive[" << id.first << " " << id.second << "]; "
+//                    << "Bbox: [" << get(m_bpm, id) << "] "
+//                    << "Point: (" << get(m_rppm, id) << ") "
+//                    << "Datum: [" << get(m_bpm, id) << "]" << std::endl;
+
+      Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+      this->insert(p);
+      return;
+    }
+
+    K2AK k2ak;
+    std::queue<AT> to_treat;
+
+    const APoint_3 ap0 = k2ak(p0);
+    const APoint_3 ap1 = k2ak(p1);
+    const APoint_3 ap2 = k2ak(p2);
+    const AFT sq_l0 = CGAL::squared_distance(ap1, ap2);
+    const AFT sq_l1 = CGAL::squared_distance(ap2, ap0);
+    const AFT sq_l2 = CGAL::squared_distance(ap0, ap1);
+
+    to_treat.push(CGAL::make_array(std::make_pair(ap0, sq_l0.sup()),
+                                   std::make_pair(ap1, sq_l1.sup()),
+                                   std::make_pair(ap2, sq_l2.sup())));
+
+    while(!to_treat.empty())
+    {
+      const AT t = std::move(to_treat.front());
+      to_treat.pop();
+
+      const APL& apl0 = t[0];
+      const APL& apl1 = t[1];
+      const APL& apl2 = t[2];
+
+      int i = (apl0.second >= apl1.second)
+                ? (apl0.second >= apl2.second) ? 0 : 2
+                : (apl1.second >= apl2.second) ? 1 : 2;
+
+      const NT max_sql = t[i].second;
+
+      // If the face is too big, do a split (two small bboxes rather than a big one)
+      if(max_sql > m_sq_length)
+      {
+        // Could be factorized, but this is simpler to read
+        if(i == 0)
+        {
+          // 0 1 2 into 0 1 A and 0 A 2
+          const APoint_3 amp = CGAL::midpoint(apl1.first, apl2.first);
+          const NT sq_half_length = apl0.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl0.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, sq_half_length),
+                                         std::make_pair(apl1.first, sq_diag_length),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, sq_diag_length),
+                                         std::make_pair(apl0.first, sq_half_length),
+                                         std::make_pair(amp, apl1.second)));
+        }
+        else if(i == 1)
+        {
+          // 0 1 2 into 0 1 A and 1 2 A
+          const APoint_3 amp = CGAL::midpoint(apl2.first, apl0.first);
+          const NT sq_half_length = apl1.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl1.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, sq_diag_length),
+                                         std::make_pair(apl1.first, sq_half_length),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, sq_half_length),
+                                         std::make_pair(apl2.first, sq_diag_length),
+                                         std::make_pair(amp, apl0.second)));
+        }
+        else // i == 2
+        {
+          // 0 1 2 into 0 A 2 and 2 A 1
+          const APoint_3 amp = CGAL::midpoint(apl0.first, apl1.first);
+          const NT sq_half_length = apl2.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl2.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, sq_half_length),
+                                         std::make_pair(apl0.first, sq_diag_length),
+                                         std::make_pair(amp, apl1.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, sq_diag_length),
+                                         std::make_pair(apl2.first, sq_half_length),
+                                         std::make_pair(amp, apl0.second)));
+        }
+      }
+      else // all edges have length below the threshold, create a primitive
+      {
+        const std::size_t pid = this->size();
+        ID id = std::make_pair(pid, fid);
+
+        put(m_dppmb, pid, data_size);
+
+        // this is basically to_double() of an APoint_3, but FT is not necessarily 'double'
+        const APoint_3& apt = compute_reference_point(apl0.first, apl1.first, apl2.first);
+        put(m_rppm, id, Point((FT(apt.x().sup()) + FT(apt.x().inf())) / FT(2),
+                              (FT(apt.y().sup()) + FT(apt.y().inf())) / FT(2),
+                              (FT(apt.z().sup()) + FT(apt.z().inf())) / FT(2)));
+
+        put(m_bpm, id, compute_bbox(apl0.first, apl1.first, apl2.first));
+
+//        std::cout << "Primitive[" << std::get<0>(id) << " " << std::get<1>(id) << " " << std::get<2>(id) << "]; "
+//                  << "Bbox: [" << get(m_bpm, id) << "] "
+//                  << "Point: (" << get(m_rppm, id) << ") "
+//                  << "Datum: [" << get(m_dpm, id) << "]" << std::endl;
+
+        Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+        this->insert(p);
+      }
+    }
+
+    ++fid;
+  }
+};
+
+} // namespace internal
+} // namespace AABB_trees
+} // namespace CGAL
+
+#endif // CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H
diff --git a/AABB_tree/test/AABB_tree/aabb_triangle_datum_covering.cpp b/AABB_tree/test/AABB_tree/aabb_triangle_datum_covering.cpp
new file mode 100644
index 00000000000..bfd756fd0ae
--- /dev/null
+++ b/AABB_tree/test/AABB_tree/aabb_triangle_datum_covering.cpp
@@ -0,0 +1,262 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/AABB_tree/internal/triangle_datum_covering.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Real_timer.h>
+
+#include <fstream>
+#include <iostream>
+
+constexpr int N = 1000;
+
+template <typename Mesh>
+void test_no_cover(const Mesh& mesh)
+{
+  std::cout << "Test WITHOUT cover" << std::endl;
+
+  using GT = typename CGAL::GetGeomTraits<Mesh>::type;
+  using FT = typename GT::FT;
+  using Point = typename GT::Point_3;
+  using Ray_3 = typename GT::Ray_3;
+  using Line_3 = typename GT::Line_3;
+
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh, CGAL::Default, CGAL::Tag_true, CGAL::Tag_true>;
+  using Traits = CGAL::AABB_traits<GT, Primitive>;
+  using Tree = CGAL::AABB_tree<Traits>;
+
+  // Build
+  Tree tree(faces(mesh).first, faces(mesh).second, mesh);
+  std::cout << faces(mesh).size() << " input faces and " << tree.size() << " primitives" << std::endl;
+
+  // Usage
+  std::cout << "Traverse..." << std::endl;
+
+  CGAL::Random r; // = CGAL::Random(1337);
+  CGAL::Random_points_in_cube_3<Point, CGAL::Creator_uniform_3<FT, Point> > g(1., r);
+
+  std::vector<Point> points;
+  points.reserve(N);
+  std::copy_n(g, N, std::back_inserter(points));
+
+  CGAL::Real_timer timer;
+  timer.start();
+
+  for(int i=0; i<N; ++i)
+  {
+    bool does_intersect = tree.do_intersect(Line_3(CGAL::ORIGIN, points[i]));
+    assert(does_intersect);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto inter_res = tree.any_intersection(Ray_3(CGAL::ORIGIN, points[i]));
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto inter_res = tree.first_intersection(Ray_3(CGAL::ORIGIN, points[i]));
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Point res = tree.closest_point(points[i]);
+    CGAL_USE(res);
+  }
+
+  for(int i=0; i<N; ++i)
+    tree.all_intersected_primitives(Line_3(CGAL::ORIGIN, points[i]), CGAL::Emptyset_iterator());
+
+  timer.stop();
+  std::cout << "Elapsed: " << timer.time() << " s." << std::endl;
+}
+
+template <typename Mesh, typename NamedParameters>
+void test_cover(const Mesh& mesh,
+                const NamedParameters& np)
+{
+  std::cout << "Test WITH cover" << std::endl;
+
+  constexpr bool check_with_standard_tree = true;
+
+  using VPM = typename CGAL::GetVertexPointMap<Mesh, NamedParameters>::const_type;
+  using GT = typename CGAL::GetGeomTraits<Mesh, NamedParameters>::type;
+  using FT = typename GT::FT;
+  using Point = typename boost::property_traits<VPM>::value_type;
+  using Ray_3 = typename GT::Ray_3;
+  using Line_3 = typename GT::Line_3;
+  using Triangle_3 = typename GT::Triangle_3;
+
+  using AABB_tree = CGAL::AABB_trees::internal::AABB_covered_triangle_tree<GT, Point>;
+
+  using FG_Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh, CGAL::Default, CGAL::Tag_true, CGAL::Tag_true>;
+  using FG_Traits = CGAL::AABB_traits<GT, FG_Primitive>;
+  using FG_Tree = CGAL::AABB_tree<FG_Traits>;
+
+  CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+
+  // Build -----------------------------------------------------------------------------------------
+  using CGAL::parameters::get_parameter;
+  using CGAL::parameters::choose_parameter;
+
+  VPM vpm = choose_parameter(get_parameter(np, CGAL::internal_np::vertex_point),
+                             get_const_property_map(CGAL::vertex_point, mesh));
+  GT gt = choose_parameter<GT>(get_parameter(np, CGAL::internal_np::geom_traits));
+
+  AABB_tree tree(diag_length / 100.);
+  tree.reserve(num_faces(mesh));
+
+  for(auto f : faces(mesh))
+  {
+    if(CGAL::Polygon_mesh_processing::is_degenerate_triangle_face(f, mesh, np))
+      continue;
+
+    const Point& p0 = get(vpm, target(halfedge(f, mesh), mesh));
+    const Point& p1 = get(vpm, target(next(halfedge(f, mesh), mesh), mesh));
+    const Point& p2 = get(vpm, source(halfedge(f, mesh), mesh));
+
+    const Triangle_3 tr = gt.construct_triangle_3_object()(p0, p1, p2);
+    tree.split_and_insert(tr);
+  }
+
+  std::cout << faces(mesh).size() << " input faces and " << tree.size() << " primitives" << std::endl;
+
+  FG_Tree fg_tree(faces(mesh).first, faces(mesh).second, mesh);
+
+  // Usage -----------------------------------------------------------------------------------------
+  std::cout << "Traverse..." << std::endl;
+
+  CGAL::Random r; // = CGAL::Random(1337);
+  CGAL::Random_points_in_cube_3<Point, CGAL::Creator_uniform_3<FT, Point> > g(1., r);
+
+  std::vector<Point> points;
+  points.reserve(N);
+  std::copy_n(g, N, std::back_inserter(points));
+
+  using AABB_traits = typename AABB_tree::AABB_traits;
+  using AABB_traversal_traits = CGAL::AABB_trees::internal::Covered_tree_traversal_traits<AABB_traits>;
+
+  using Do_intersect_traits = typename AABB_traversal_traits::template Do_intersect_traits<Line_3>;
+  using First_intersection_traits = typename AABB_traversal_traits::template First_intersection_traits<Ray_3>;
+  using Projection_traits = typename AABB_traversal_traits::Projection_traits;
+
+  using size_type = typename AABB_tree::size_type;
+  using Primitive_id = typename AABB_tree::Primitive_id;
+  using Counting_iterator = CGAL::internal::AABB_tree::Counting_output_iterator<Primitive_id, size_type>;
+  using Listing_primitive_traits = typename AABB_traversal_traits::template Listing_primitive_traits<Line_3, Counting_iterator>;
+
+  CGAL::Real_timer timer;
+  timer.start();
+
+  for(int i=0; i<N; ++i)
+  {
+    Do_intersect_traits do_intersect_traversal_traits(tree.traits());
+    Line_3 query(CGAL::ORIGIN, points[i]);
+    tree.traversal(query, do_intersect_traversal_traits);
+    bool does_intersect = do_intersect_traversal_traits.is_intersection_found();
+    assert(does_intersect);
+
+    if(check_with_standard_tree)
+      assert(does_intersect == fg_tree.do_intersect(query));
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Ray_3 query(CGAL::ORIGIN, points[i]);
+    First_intersection_traits first_intersection_traversal_traits(tree.traits());
+    tree.traversal(query, first_intersection_traversal_traits);
+    auto inter_res = first_intersection_traversal_traits.result();
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Ray_3 query(CGAL::ORIGIN, points[i]);
+    auto inter_res = tree.first_intersection(query);
+    assert(inter_res);
+
+    if(check_with_standard_tree)
+    {
+      auto fg_inter_res = fg_tree.first_intersection(query);
+      assert(inter_res->first == fg_inter_res->first);
+    }
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto hint = tree.best_hint(points[i]);
+    Projection_traits projection_traits(hint.first, hint.second, tree.traits());
+    Point res = tree.closest_point(points[i]);
+
+    if(check_with_standard_tree)
+    {
+      Point fg_res = fg_tree.closest_point(points[i]);
+      assert(res == fg_res);
+    }
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Line_3 query(CGAL::ORIGIN, points[i]);
+
+    size_type counter = 0;
+    Counting_iterator out(&counter);
+    Listing_primitive_traits listing_traversal_traits(out, tree.traits());
+    tree.traversal(query, listing_traversal_traits);
+    assert(counter > 0);
+
+    if(check_with_standard_tree)
+    {
+      typename FG_Tree::size_type fg_counter = 0;
+      CGAL::internal::AABB_tree::Counting_output_iterator<typename FG_Tree::Primitive_id,
+                                                          typename FG_Tree::size_type> fg_out(&fg_counter);
+      fg_tree.all_intersected_primitives(query, fg_out);
+      assert(counter == fg_counter);
+    }
+  }
+
+  timer.stop();
+  std::cout << "Elapsed: " << timer.time() << " s." << std::endl;
+}
+
+template <typename K>
+void test(const char* filename)
+{
+  std::cout << " === Test Kernel " << typeid(K).name() << " === " << std::endl;
+
+  using Point_3 = typename K::Point_3;
+  using Mesh = CGAL::Surface_mesh<Point_3>;
+
+  Mesh mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh) ||
+     !is_triangle_mesh(mesh))
+  {
+    std::cerr << "Error: invalid input" << std::endl;
+    return;
+  }
+
+  test_no_cover(mesh);
+  test_cover(mesh, CGAL::parameters::default_values());
+}
+
+int main(int, char**)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  test<CGAL::Exact_predicates_inexact_constructions_kernel>("data/bunny00.off");
+  test<CGAL::Exact_predicates_exact_constructions_kernel>("data/bunny00.off");
+}
+
diff --git a/Intersections_2/include/CGAL/Intersection_traits.h b/Intersections_2/include/CGAL/Intersection_traits.h
index 137bb500035..96c28e8e899 100644
--- a/Intersections_2/include/CGAL/Intersection_traits.h
+++ b/Intersections_2/include/CGAL/Intersection_traits.h
@@ -88,28 +88,6 @@ struct Intersection_traits {
 };
 
 
-// Alias that gets the Kernel automatically and does some error checking.
-// Including corresponding specialization for Bbox, as it has no Kernel.
-template<typename A, typename B>
-class IT : public Intersection_traits< typename Kernel_traits<A>::Kernel, A, B > {
-  typedef typename Kernel_traits<A>::Kernel A_Kernel;
-  typedef typename Kernel_traits<B>::Kernel B_Kernel;
-  // CGAL_static_assertion_msg( (boost::is_same< A_Kernel, B_Kernel>::value),
-  //                            "IT instantiated with objects from two different Kernels");
-};
-
-class Bbox_2;
-class Bbox_3;
-
-template<typename B>
-class IT<Bbox_2, B> : public Intersection_traits< typename Kernel_traits<B>::Kernel, CGAL::Bbox_2, B >
-{ };
-
-template<typename B>
-class IT<Bbox_3, B> : public Intersection_traits< typename Kernel_traits<B>::Kernel, CGAL::Bbox_3, B >
-{ };
-
-
 namespace Intersections {
 namespace internal {
 
diff --git a/Intersections_2/include/CGAL/Intersection_traits_2.h b/Intersections_2/include/CGAL/Intersection_traits_2.h
index d87c599548c..921405c27f9 100644
--- a/Intersections_2/include/CGAL/Intersection_traits_2.h
+++ b/Intersections_2/include/CGAL/Intersection_traits_2.h
@@ -14,6 +14,7 @@
 #define CGAL_INTERSECTION_TRAITS_2_H
 
 #include <CGAL/Intersection_traits.h>
+#include <CGAL/Bbox_2.h>
 
 #include <boost/variant.hpp>
 #include <boost/optional.hpp>
diff --git a/Intersections_3/include/CGAL/Intersection_traits_3.h b/Intersections_3/include/CGAL/Intersection_traits_3.h
index 84cac64e4bd..85287e8c044 100644
--- a/Intersections_3/include/CGAL/Intersection_traits_3.h
+++ b/Intersections_3/include/CGAL/Intersection_traits_3.h
@@ -20,6 +20,7 @@
 #define CGAL_INTERSECTION_TRAITS_3_H
 
 #include <CGAL/Intersection_traits.h>
+#include <CGAL/Bbox_3.h>
 
 #include <boost/optional.hpp>
 #include <boost/variant.hpp>