diff --git a/Convex_hull_3/benchmark/Convex_hull_3/CMakeLists.txt b/Convex_hull_3/benchmark/Convex_hull_3/CMakeLists.txt
new file mode 100644
index 000000000000..6c407de138b4
--- /dev/null
+++ b/Convex_hull_3/benchmark/Convex_hull_3/CMakeLists.txt
@@ -0,0 +1,18 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.1...3.22)
+project(Bench_convex_hull_3)
+
+find_package(CGAL REQUIRED)
+
+include_directories(BEFORE "include")
+
+# 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()
\ No newline at end of file
diff --git a/Convex_hull_3/benchmark/Convex_hull_3/bench_ch3_do_intersect.cpp b/Convex_hull_3/benchmark/Convex_hull_3/bench_ch3_do_intersect.cpp
new file mode 100644
index 000000000000..1a0eeda8b229
--- /dev/null
+++ b/Convex_hull_3/benchmark/Convex_hull_3/bench_ch3_do_intersect.cpp
@@ -0,0 +1,356 @@
+#define CGAL_PROFILE_CONVEX_HULL_DO_INTERSECT
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Extreme_points_traits_adapter_3.h>
+#include <CGAL/convex_hull_3.h>
+#include <CGAL/Convex_hull_hierarchy_3.h>
+#include <CGAL/Convex_hull_3/do_intersect.h>
+#include <CGAL/Convex_hull_3/distance.h>
+
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <CGAL/optimal_bounding_box.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
+
+#include <vector>
+#include <fstream>
+
+template <typename K>
+struct Test
+{
+  typedef typename K::FT              FT;
+  typedef typename K::Point_3         P;
+  typedef typename K::Segment_3       S;
+  typedef typename K::Vector_3        V;
+  typedef typename K::Ray_3           R;
+  typedef typename K::Line_3          L;
+  typedef typename K::Triangle_3      T;
+  typedef typename K::Plane_3         Pl;
+  typedef typename K::Tetrahedron_3   Tet;
+  typedef typename K::Iso_cuboid_3    Cub;
+  typedef typename std::vector<P>     PR;
+
+private:
+  CGAL::Random& r;
+  double m = 0, M = 1;
+  const double PI=3.14159265358979323846;
+
+public:
+   Test(CGAL::Random& r) : r(r) { }
+
+private:
+  P random_point() const
+  {
+    return P(FT(r.get_double(m, M)), FT(r.get_double(m, M)), FT(r.get_double(m, M)));
+  }
+
+  P random_sphere_point() const
+  {
+    FT a=r.get_double(0,2*PI);
+    FT b=r.get_double(-PI/2,PI/2);
+    return P( std::cos(a)*std::cos(b), std::sin(a)*std::cos(b), std::sin(b));
+  }
+
+  void test_Tet(std::vector< PR > &a, std::vector< PR > &b)
+  {
+    std::cout << "Do intersect of package Intersection_3" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(std::size_t i=0; i<a.size(); ++i)
+      CGAL::do_intersect(Tet(a[i][0], a[i][1], a[i][2], a[i][3]), Tet(b[i][0], b[i][1], b[i][2], b[i][3]));
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+    std::cout << "Do intersect of Convex_hull" << std::endl;
+    t.reset();
+    t.start();
+    for(std::size_t i=0; i<a.size(); ++i)
+      CGAL::Convex_hull_3::do_intersect(a[i], b[i]);
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+    std::cout << "Separation distance of Convex_hull" << std::endl;
+    t.reset();
+    t.start();
+    for(std::size_t i=0; i<a.size(); ++i)
+      CGAL::Convex_hull_3::experimental::separation_distance(a[i], b[i]);
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl << std::endl;
+  }
+
+  void Tet_tet(int N)
+  {
+    std::cout << "Tetrahedron likely to intersect" << std::endl;
+    std::vector<PR> a;
+    std::vector<PR> b;
+    for(int i=0; i<N; ++i)
+      a.push_back(PR({random_point(), random_point(), random_point(), random_point()}));
+    for(int i=0; i<N; ++i)
+      b.push_back(PR({random_point(), random_point(), random_point(), random_point()}));
+    test_Tet(a, b);
+  }
+
+  void Tet_mirror(int N)
+  {
+    std::cout << "Tetrahedrons closed" << std::endl;
+    std::vector<PR> a;
+    std::vector<PR> b;
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      if(p0.x() > p0.y())
+        p0 = P(-p0.x(), -p0.y(), -p0.z());
+      P p1 = random_point();
+      if(p1.x() > p1.y())
+        p1 = P(-p1.x(), -p1.y(), -p1.z());
+      P p2 = random_point();
+      if(p2.x() > p2.y())
+        p2 = P(-p2.x(), -p2.y(), -p2.z());
+      P p3 = random_point();
+      if(p3.x() > p3.y())
+        p3 = P(-p3.x(), -p3.y(), -p3.z());
+
+      P q0 = random_point();
+      if(q0.x() < q0.y())
+        q0 = P(-q0.x(), -q0.y(), -q0.z());
+      P q1 = random_point();
+      if(q1.x() < q1.y())
+        q1 = P(-q1.x(), -q1.y(), -q1.z());
+      P q2 = random_point();
+      if(q2.x() < q2.y())
+        q2 = P(-q2.x(), -q2.y(), -q2.z());
+      P q3 = random_point();
+      if(q3.x() < q3.y())
+        q3 = P(-q3.x(), -q3.y(), -q3.z());
+
+      a.push_back(PR({p0,p1,p2,p3}));
+      b.push_back(PR({q0,q1,q2,q3}));
+    }
+    test_Tet(a, b);
+  }
+
+  void Tet_stretched(int N)
+  {
+    std::cout << "Tetrahedrons streched" << std::endl;
+    std::vector<PR> a;
+    std::vector<PR> b;
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point()+V(-1,1,0);
+      P p1 = random_point()+V(-1,1,0);
+      P p2 = random_point()+V(1,-1,0);
+      P p3 = random_point()+V(1,-1,0);
+
+      P q0 = random_point()+V(-1,-1,0);
+      P q1 = random_point()+V(-1,-1,0);
+      P q2 = random_point()+V(1,1,0);
+      P q3 = random_point()+V(1,1,0);
+
+      a.push_back(PR({p0,p1,p2,p3}));
+      b.push_back(PR({q0,q1,q2,q3}));
+    }
+    test_Tet(a, b);
+  }
+
+  void Tet_shift(int N)
+  {
+    std::cout << "Tetrahedrons that do not intersect" << std::endl;
+    std::vector<PR> a;
+    std::vector<PR> b;
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P p3 = random_point();
+
+      P q0 = random_point()+V(1,1,1);
+      P q1 = random_point()+V(1,1,1);
+      P q2 = random_point()+V(1,1,1);
+      P q3 = random_point()+V(1,1,1);
+
+      a.push_back(PR({p0,p1,p2,p3}));
+      b.push_back(PR({q0,q1,q2,q3}));
+    }
+    test_Tet(a, b);
+
+  }
+
+  void Test_sphere(int N, int M)
+  {
+    std::cout << "Sphere of size " << M << std::endl;
+    CGAL::Real_timer t;
+    std::vector< std::vector<P> > points_ranges;
+    std::vector< CGAL::Surface_mesh<P> > hulls;
+    std::vector< CGAL::Convex_hull_hierarchy_3< CGAL::Surface_mesh<P> > > hulls_hierarchy;
+
+    double nb_tests=double(N*(N-1)/2);
+
+    for(int i=0; i<N; ++i)
+    {
+      std::vector<P> sp;
+      V vec( P(0,0,0), random_sphere_point());
+      vec*=1.2;
+      for(int i=0; i<M; ++i)
+        sp.push_back(random_sphere_point()+vec);
+      points_ranges.push_back(std::move(sp));
+    }
+
+    std::cout << "Compute convex hulls" << std::endl;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      CGAL::Surface_mesh<P> hull;
+      CGAL::convex_hull_3(points_ranges[i].begin(), points_ranges[i].end(), hull);
+      hulls.push_back(std::move(hull));
+    }
+    t.stop();
+    std::cout << "  " << t.time()/N << " sec" << std::endl;
+    t.reset();
+
+    std::cout << "Construct hierarchy of convex hulls" << std::endl;
+    t.start();
+    for(int i=0; i<N; ++i)
+      hulls_hierarchy.emplace_back(hulls[i], CGAL::parameters::compute_convex_hull(false));
+    t.stop();
+    double time_build = t.time()/N;
+    std::cout << "  " << t.time()/N << " sec" << std::endl;
+    t.reset();
+
+    std::cout << "Do intersect with RangePoints" << std::endl;
+    nb_visited=0;
+    t.start();
+    for(int i=0; i<N; ++i)
+      for(int j=i+1; j<N; ++j)
+        CGAL::Convex_hull_3::do_intersect(points_ranges[i], points_ranges[j]);
+    t.stop();
+
+    std::cout << "  " << t.time()/nb_tests << " sec per test" << std::endl;
+    t.reset();
+
+    std::cout << "Do intersect with Hulls" << std::endl;
+    nb_visited=0;
+    t.start();
+    for(int i=0; i<N; ++i)
+      for(int j=i+1; j<N; ++j)
+        CGAL::Convex_hull_3::do_intersect(hulls[i], hulls[j]);
+    t.stop();
+    double time_without_hierarchy = t.time()/nb_tests;
+    std::cout << "  " << t.time()/nb_tests << " sec per test" << std::endl;
+    t.reset();
+
+    std::cout << "Do intersect with Hierarchical Hulls" << std::endl;
+    nb_visited=0;
+    t.start();
+    for(int i=0; i<N; ++i)
+      for(int j=i+1; j<N; ++j)
+        CGAL::Convex_hull_3::do_intersect(hulls_hierarchy[i], hulls_hierarchy[j]);
+    t.stop();
+
+    std::cout << "  " << t.time()/nb_tests << " sec per test" << std::endl;
+    std::cout << "     " << time_build / (time_without_hierarchy - t.time()/nb_tests) << " ratio tests/build";
+    t.reset();
+    std::cout << std::endl;
+  }
+
+public:
+  void run()
+  {
+    std::cout << "Kernel: " << typeid(K).name() << std::endl;
+    Tet_tet(1000);
+    Tet_mirror(1000);
+    Tet_shift(1000);
+    Tet_stretched(1000);
+
+    Test_sphere(1000, 40);
+    Test_sphere(1000, 100);
+    Test_sphere(1000, 350);
+    Test_sphere(1000, 1000);
+    Test_sphere(300, 3500);
+    Test_sphere(50, 10000);
+    Test_sphere(25, 35000);
+    Test_sphere(25, 100000);
+    Test_sphere(10, 350000);
+    Test_sphere(10, 1000000);
+    std::cout << std::endl;
+  }
+};
+
+template <typename K>
+void bench_on_data(const std::string &f1, const std::string &f2){
+  typedef typename K::Point_3                                               Point_3;
+  typedef CGAL::Surface_mesh<Point_3>                              Mesh;
+
+  std::vector<typename K::Point_3> pts1, pts2;
+  std::vector<boost::container::small_vector<std::size_t, 3>> trs1, trs2;
+  if (!CGAL::IO::read_polygon_soup(f1, pts1, trs1))
+  {
+    std::cerr << "Cannot read " << f1 << "\n";
+  }
+  if (!CGAL::IO::read_polygon_soup(f2, pts2, trs2))
+  {
+    std::cerr << "Cannot read " << f2 << "\n";
+  }
+
+  CGAL::Real_timer t;
+  t.start();
+  Mesh hull1, hull2;
+  CGAL::convex_hull_3(pts1.begin(), pts1.end(), hull1);
+  CGAL::convex_hull_3(pts2.begin(), pts2.end(), hull2);
+  t.stop();
+  std::cout << "Computing convex hulls: " << t.time() << " sec" << std::endl;
+  std::cout << "Convex hull size:" << vertices(hull1).size() << ", " << vertices(hull2).size() << "\n" << std::endl;
+  t.reset();
+
+  t.start();
+  CGAL::Convex_hull_3::do_intersect(hull1, hull2);
+  t.stop();
+  std::cout << "Do intersect: " << t.time() << " sec\n" << std::endl;
+  t.reset();
+
+  t.start();
+  std::array<Point_3, 8> obb1, obb2;
+  CGAL::oriented_bounding_box(hull1, obb1, CGAL::parameters::use_convex_hull(false));
+  CGAL::oriented_bounding_box(hull2, obb2, CGAL::parameters::use_convex_hull(false));
+  t.stop();
+  std::cout << "Computing Obbs: " << t.time() << " sec\n" << std::endl;
+  t.reset();
+
+  t.start();
+  CGAL::Convex_hull_3::do_intersect(obb1, obb2);
+  t.stop();
+  std::cout << "Do intersect with Obbs: " << t.time() << " sec\n" << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+  // std::cout.precision(17);
+  // std::cerr.precision(17);
+
+  std::cout << "3D Distance tests" << std::endl;
+
+  CGAL::Random rp;
+  CGAL::Random r(argc==1?rp.get_seed():std::stoi(argv[1]));
+  // std::cout << "random seed = " << r.get_seed() << std::endl;
+
+  // Test<CGAL::Simple_cartesian<double> >(r).run();
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r).run();
+  // Test<CGAL::Exact_predicates_exact_constructions_kernel>(r).run();
+
+  const std::string f1 = (argc>2) ? argv[2] : CGAL::data_file_path("meshes/elephant.off");
+  const std::string f2 = (argc>3) ? argv[3] : CGAL::data_file_path("meshes/sphere.off");
+  // bench_on_data<CGAL::Exact_predicates_inexact_constructions_kernel>(f1,f2);
+
+  std::cout << "Done!" << std::endl;
+
+  return EXIT_SUCCESS;
+}
diff --git a/Convex_hull_3/doc/Convex_hull_3/CGAL/convex_hull_3.h b/Convex_hull_3/doc/Convex_hull_3/CGAL/convex_hull_3.h
index f6621de0063a..d6bbd232db7e 100644
--- a/Convex_hull_3/doc/Convex_hull_3/CGAL/convex_hull_3.h
+++ b/Convex_hull_3/doc/Convex_hull_3/CGAL/convex_hull_3.h
@@ -154,9 +154,9 @@ then the default traits class used is `Convex_hull_traits_3<R>`, and `R` otherwi
 */
 template <class InputRange, class OutputIterator, class Traits>
 OutputIterator
-extreme_points_3(InputRange range,
-              OutputIterator out,
-              const Traits& traits);
+extreme_points_3(const InputRange &range,
+                 OutputIterator out,
+                 const Traits& traits);
 
 
 
diff --git a/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt b/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt
index 0722e7456fb8..51a890f55fcd 100644
--- a/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt
+++ b/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt
@@ -6,7 +6,7 @@ namespace CGAL {
 \anchor Chapter_3D_Convex_Hulls
 \anchor chapconvexhull3
 \cgalAutoToc
-\authors Susan Hert and Stefan Schirra
+\authors Samuel Hornus, Susan Hert, Stefan Schirra, and Leo Valque
 
 \cgalFigureBegin{figchbimba,chull_bimba.png}
 The convex hull of a model made of 192135 points.
@@ -37,7 +37,7 @@ The function
 implementation of the quickhull algorithm \cgalCite{bdh-qach-96}.
 There are two versions of this
 function available, one that can be used when it is known that the output
-will be a polyhedron (<I>i.e.</I>, there are more than three points and
+will be a polyhedron (i.e., there are more than three points and
 they are not all collinear) and one that handles all degenerate cases
 and returns an `Object`, which may be a point, a segment, a
 triangle, or a polyhedron. Both versions accept a range of input
@@ -84,7 +84,7 @@ Note that the latter may also be planar polygon with a border.
 
 \cgalExample{Convex_hull_3/quickhull_any_dim_3.cpp}
 
-\subsection Convex_hull_3ExtremePoints Extreme points
+\subsection Convex_hull_3ExtremePoints Extreme Points
 In addition to the `convex_hull_3()` function, the function `extreme_points_3()` is also
 provided in case only the points on the convex hull are required (without the connectivity
 information). In addition the traits class adapter `CGAL::Extreme_points_traits_adapter_3`
@@ -133,7 +133,7 @@ point set or not. This function is used in postcondition testing for
 
 Fully dynamic maintenance of a convex hull can be achieved by using the
 class `Delaunay_triangulation_3`. This class supports insertion
-and removal of points (<I>i.e.</I>, vertices of the triangulation) and the
+and removal of points (i.e., vertices of the triangulation) and the
 convex hull edges are simply the finite edges of infinite faces.
 The following example illustrates the dynamic construction of a convex hull.
 First, random points from a sphere of a certain radius are generated and are
@@ -154,6 +154,18 @@ that is model of the concept `MutableFaceGraph`, e.g. `Polyhedron_3` and `Surfac
 
 \cgalExample{Convex_hull_3/dynamic_hull_3.cpp}
 
+\section Convex_hull_3Do_intersect Intersection Test of Convex Hulls
+
+`CGAL::Convex_hull_3::do_intersect(ch1, ch2)` determines whether two convex hulls intersect or not.
+
+Input hulls can be provided as a range of points or as a graph, and may be of different types.
+Furthermore, when many intersection queries use the same object, one should wrap that input in the class `CGAL::Convex_hull_hierarchy_3` as to construct an optimized view of the convex hull and accelerate intersection tests.
+This is especially true when the convex hull is made of a large number of vertices (see `CGAL::Convex_hull_hierarchy_3` for more details).
+
+\subsection Do_intersectExample Example
+The following program reads points from two input files, computes their convex hulls with `Convex_hull_hierarchy_3` and prints if they intersect or not.
+\cgalExample{Convex_hull_3/do_intersect_ch3.cpp}
+
 \section Convex_hull_3Performance Performance
 
 In the following, we compare the running times of the two approaches to compute 3D convex hulls.
@@ -169,6 +181,13 @@ the static approach needed 0.18s, while the dynamic approach needed 1.90s.
 The measurements have been performed using \cgal 3.9, using the \gnu \cpp compiler version 4.3.5, under Linux (Debian distribution),
 with the compilation options <TT>-O3 -DCGAL_NDEBUG</TT>. The computer used was equipped with a 64bit Intel Xeon 2.27GHz processor and 12GB of RAM.
 
+
+\section Convex_hull_3History Implementation History
+
+The initial version of this package was implemented by Susan Hert and Stefan Schirra, and released in 2011.
+The functions `do_intersect()`, `extreme_point_3()`, and the class `Convex_hull_hierarchy_3` were implemented in 2026 by Samuel Hornus and Leo Valque, under the supervision of Sébastien Loriot.
+
+
 */
 } /* namespace CGAL */
 
diff --git a/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in b/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in
index 83e604ed3ba4..1056c8343925 100644
--- a/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in
+++ b/Convex_hull_3/doc/Convex_hull_3/Doxyfile.in
@@ -1,4 +1,11 @@
 @INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
 
+
+EXCLUDE_SYMBOLS += experimental
+
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Convex Hulls"
 INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Convex_hull_3/distance.h
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Convex_hull_3/do_intersect.h
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Convex_hull_hierarchy_3.h
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/extreme_point_3.h
diff --git a/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt b/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt
index f8530e90c99e..6f6ad88f3c89 100644
--- a/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt
+++ b/Convex_hull_3/doc/Convex_hull_3/PackageDescription.txt
@@ -6,9 +6,15 @@
 /// \defgroup PkgConvexHull3Traits Traits Classes
 /// \ingroup PkgConvexHull3Ref
 
+/// \defgroup PkgConvexHull3Intersections Intersection Test Functions
+/// \ingroup PkgConvexHull3Ref
+
+/// \defgroup PkgConvexHull3Queries Directional Queries Functions
+/// \ingroup PkgConvexHull3Ref
+
 /*! \defgroup PkgConvexHull3Functions Convex Hull Functions
    \ingroup PkgConvexHull3Ref
-The function `convex_hull_3()` computes the convex hull of a given set of
+The function `CGAL::convex_hull_3()` computes the convex hull of a given set of
 three-dimensional points.
 
 Two versions of this function are available. The first can be used when it
@@ -27,7 +33,7 @@ degenerate hull may also be possible.
 \cgalPkgPicture{Convex_hull_3/fig/bunny.png}
 \cgalPkgSummaryBegin
 \cgalPkgAuthors{Susan Hert and Stefan Schirra}
-\cgalPkgDesc{This package provides functions  for computing convex hulls in three dimensions as well as functions for checking if sets of points are strongly convex or not. One can compute the convex hull of a set of points in three dimensions in two ways: using a static algorithm or using a triangulation to get a fully dynamic computation.}
+\cgalPkgDesc{This package provides functions for computing convex hulls in three dimensions, as well as functions for checking if sets of points are strongly convex or not, and checking if two convex hulls intersect or not. One can compute the convex hull of a set of points in three dimensions in three ways: using a static algorithm, using a triangulation to get a fully dynamic computation, or using the specific data structure `Convex_hull_hierarchy_3` for optimal intersection testing.}
 \cgalPkgManuals{Chapter_3D_Convex_Hulls,PkgConvexHull3Ref}
 \cgalPkgSummaryEnd
 \cgalPkgShortInfoBegin
@@ -48,7 +54,7 @@ polytope with vertices in `P`.  A point in `P` is an extreme point
 
 \cgal provides functions for computing convex hulls in two, three
 and arbitrary dimensions as well as functions for testing if a given set of
-points in is strongly convex or not.  This chapter describes the functions
+points in is strongly convex or not and checking if two convex hulls intersect. This chapter describes the functions
 available for three dimensions.
 
 \cgalCRPSection{Assertions}
@@ -69,12 +75,19 @@ defining `CGAL_CHECK_EXPENSIVE`.
 
 - `CGAL::Convex_hull_traits_3<R>`
 - `CGAL::Extreme_points_traits_adapter_3`
+- `CGAL::make_extreme_points_traits_adapter()`
 
 \cgalCRPSection{Convex Hull Functions}
 
-- `CGAL::convex_hull_3`
-- `CGAL::extreme_points_3`
-- `CGAL::make_extreme_points_traits_adapter`
+- `CGAL::convex_hull_3()`
+- `CGAL::extreme_points_3()`
+
+\cgalCRPSection{Intersection and Directional Queries}
+
+- `CGAL::extreme_point_3()`
+- `CGAL::extreme_vertex_3()`
+- `CGAL::Convex_hull_3::do_intersect()`
+- `CGAL::Convex_hull_hierarchy_3`
 
 \cgalCRPSection{Convexity Checking Function}
 
diff --git a/Convex_hull_3/doc/Convex_hull_3/dependencies b/Convex_hull_3/doc/Convex_hull_3/dependencies
index d4ca856283cf..7f0207c28d2a 100644
--- a/Convex_hull_3/doc/Convex_hull_3/dependencies
+++ b/Convex_hull_3/doc/Convex_hull_3/dependencies
@@ -7,6 +7,7 @@ Stream_support
 Convex_hull_2
 Convex_hull_d
 Triangulation_3
+Property_map
 Polyhedron
 Surface_mesh
 BGL
diff --git a/Convex_hull_3/doc/Convex_hull_3/examples.txt b/Convex_hull_3/doc/Convex_hull_3/examples.txt
index ed7a90431515..0cab7cdfddeb 100644
--- a/Convex_hull_3/doc/Convex_hull_3/examples.txt
+++ b/Convex_hull_3/doc/Convex_hull_3/examples.txt
@@ -8,4 +8,5 @@
 \example Convex_hull_3/extreme_points_3_sm.cpp
 \example Convex_hull_3/extreme_indices_3.cpp
 \example Convex_hull_3/quickhull_indexed_triangle_set_3.cpp
+\example Convex_hull_3/do_intersect_ch3.cpp
 */
diff --git a/Convex_hull_3/examples/Convex_hull_3/CMakeLists.txt b/Convex_hull_3/examples/Convex_hull_3/CMakeLists.txt
index 99b49365013e..5f6a30688ac3 100644
--- a/Convex_hull_3/examples/Convex_hull_3/CMakeLists.txt
+++ b/Convex_hull_3/examples/Convex_hull_3/CMakeLists.txt
@@ -18,6 +18,7 @@ create_single_source_cgal_program("graph_hull_3.cpp")
 create_single_source_cgal_program("quickhull_any_dim_3.cpp")
 create_single_source_cgal_program("extreme_points_3_sm.cpp")
 create_single_source_cgal_program("extreme_indices_3.cpp")
+create_single_source_cgal_program("do_intersect_ch3.cpp")
 
 find_package(OpenMesh QUIET)
 if(OpenMesh_FOUND)
diff --git a/Convex_hull_3/examples/Convex_hull_3/do_intersect_ch3.cpp b/Convex_hull_3/examples/Convex_hull_3/do_intersect_ch3.cpp
new file mode 100644
index 000000000000..c7cc41599312
--- /dev/null
+++ b/Convex_hull_3/examples/Convex_hull_3/do_intersect_ch3.cpp
@@ -0,0 +1,51 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Extreme_points_traits_adapter_3.h>
+#include <CGAL/convex_hull_3.h>
+#include <CGAL/Convex_hull_3/distance.h>
+#include <CGAL/Convex_hull_3/do_intersect.h>
+#include <CGAL/Convex_hull_hierarchy_3.h>
+
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <vector>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel      K;
+typedef K::Point_3                                               Point_3;
+typedef CGAL::Surface_mesh<Point_3>                              Mesh;
+typedef Mesh::Property_map<Mesh::Vertex_index, Point_3>          PointMap;
+typedef CGAL::Convex_hull_hierarchy_3<Mesh>                      Convex_hull_hierarchy_3;
+
+int main(int argc, char* argv[])
+{
+  const std::string f1 = (argc>1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+  const std::string f2 = (argc>2) ? argv[2] : CGAL::data_file_path("meshes/sphere.off");
+  Mesh sm1, sm2;
+  if(!CGAL::IO::read_polygon_mesh(f1, sm1))
+  {
+    std::cerr<< "Cannot open " << f1 <<std::endl;
+    return 1;
+  }
+  if(!CGAL::IO::read_polygon_mesh(f2, sm2))
+  {
+    std::cerr<< "Cannot open " << f2 <<std::endl;
+    return 1;
+  }
+
+  //call the function with the traits adapter for vertices
+  bool res = CGAL::Convex_hull_3::do_intersect(vertices(sm1), vertices(sm2),
+                                               CGAL::parameters::point_map(sm1.points()),
+                                               CGAL::parameters::point_map(sm2.points()));
+
+  // When performing many intersection tests, we can construct an optimized version of the convex hull with the class `Convex_hull_hierarchy_3`
+  Convex_hull_hierarchy_3 hsm1(sm1);
+  Convex_hull_hierarchy_3 hsm2(sm2);
+  res = CGAL::Convex_hull_3::do_intersect(hsm1, hsm2);
+  std::cout << "do convex hulls intersect? " << std::boolalpha << res << "\n";
+
+
+  return 0;
+}
\ No newline at end of file
diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/distance.h b/Convex_hull_3/include/CGAL/Convex_hull_3/distance.h
new file mode 100644
index 000000000000..3ef38c829f37
--- /dev/null
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/distance.h
@@ -0,0 +1,430 @@
+// Copyright (c) 2022 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_CONVEX_HULL_3_DISTANCE_H
+#define CGAL_CONVEX_HULL_3_DISTANCE_H
+
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
+#include <CGAL/Kernel_23/internal/Has_boolean_tags.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Convex_hull_hierarchy_3.h>
+#include <CGAL/Convex_hull_3/internal/helpers.h>
+#include <CGAL/squared_distance_3.h>
+
+#include <CGAL/Container_helper.h>
+#include <CGAL/Dynamic_property_map.h>
+#include <CGAL/property_map.h>
+
+#include <boost/range/value_type.hpp>
+#include <boost/range/reference.hpp>
+
+#include <boost/graph/adjacency_list.hpp>
+
+#include <CGAL/Real_timer.h>
+
+#include <vector>
+
+namespace CGAL {
+namespace Convex_hull_3 {
+namespace experimental {
+
+/*
+Note: The function is not stable with EPICK (possible infinite loops), It needs to be rewritten by moving all constructions inside predicates and by defining objects (Vector_3, Direction_3,...) are defined directly from the input points.
+For instance, Vector_3 should be defined by two points and Direction_3 by 3 vectors (i.e 6 points) if the simplex size is 3 or 4.
+Currently the function simply switches to EPECK and is considered experimental.
+*/
+
+namespace predicates_impl {
+
+// returns true iff pt is on the negative side of the plane defined by (ep0, ep1) and normal
+template <class K>
+inline bool
+strictly_on_left_of_triangle_edge(const typename K::Vector_3& pt,
+                                  const typename K::Vector_3& normal,
+                                  const typename K::Vector_3& edge,
+                                  const K& k)
+{
+  const typename K::Construct_cross_product_vector_3 cross=k.construct_cross_product_vector_3_object();
+  const typename K::Compute_scalar_product_3 dot=k.compute_scalar_product_3_object();
+
+  return is_negative(dot(cross(edge, normal), pt));
+}
+
+/*
+Compute the Voronoi cell containing the origin and return the direction of the simplex closest to it.
+The origin necessarily lies below the edge formed by the first two vertices of the triangle,
+we exploit this fact to simplify the computation.
+*/
+template<typename K, typename Vector_3>
+Vector_3 triangle_dir_to_origin(boost::container::small_vector<Vector_3, 4>& simplex, const K &k){
+  using FT=typename K::FT;
+
+  const typename K::Construct_cross_product_vector_3 cross=k.construct_cross_product_vector_3_object();
+  const typename K::Compute_scalar_product_3 dot=k.compute_scalar_product_3_object();
+
+  Vector_3 &a = simplex[0];
+  Vector_3 &b = simplex[1];
+  Vector_3 &c = simplex[2];
+  Vector_3 ab = b-a;
+  Vector_3 ca = a-c;
+  Vector_3 bc = c-b;
+
+  Vector_3 n=cross(ca, bc);
+  assert(n!=NULL_VECTOR);
+
+  // FT ratio_ab=dot(ab, -a)/ab.squared_length();
+  FT ratio_ca=dot(ca, -c);
+  FT ratio_cb=dot(-bc, -c);
+
+  bool on_left_of_ca = strictly_on_left_of_triangle_edge(-c, n, ca, k);
+  bool on_left_of_bc = strictly_on_left_of_triangle_edge(-b, n, bc, k);
+
+  // The origin cannot be on the side of ab, we abuse of this
+  if(/*!on_left_of_ab &&*/ !on_left_of_ca && !on_left_of_bc){
+    if(is_negative(dot(n, -a))){
+      return -n;
+    }
+    std::swap(simplex[1], simplex[2]);
+    return n;
+  }
+
+  if(on_left_of_ca && ratio_ca>=0 /*&& ratio_ca<=1*/){
+    Vector_3 dir=-cross(cross(ca, c), ca);
+    simplex[1] = simplex[2];
+    simplex.pop_back();
+    return dir;
+  }
+  if(on_left_of_bc && ratio_cb>=0 /*&& ratio_cb<=1*/){
+    Vector_3 dir=-cross(cross(bc, b), bc);
+    simplex[0] = simplex[2];
+    simplex.pop_back();
+    return dir;
+  }
+
+  simplex[0] = simplex[2];
+  simplex.pop_back();
+  simplex.pop_back();
+  return -simplex[0];
+}
+
+/*
+Compute the Voronoi cell containing the origin and return the direction of the simplex closest to it.
+The origin necessarily lies below the triangle formed by the first three vertices of the tetrahedron,
+we exploit this fact to simplify the computation.
+*/
+template<typename K, typename Vector_3>
+Vector_3 tetrahedron_dir_to_origin(boost::container::small_vector<Vector_3, 4>& simplex, const K &k) {
+  const typename K::Construct_cross_product_vector_3 cross=k.construct_cross_product_vector_3_object();
+  const typename K::Compute_scalar_product_3 dot=k.compute_scalar_product_3_object();
+  const typename K::Orientation_3 orientation = k.orientation_3_object();
+
+  Vector_3 &a=simplex[0];
+  Vector_3 &b=simplex[1];
+  Vector_3 &c=simplex[2];
+  Vector_3 &d=simplex[3];
+
+  Vector_3 ad=d-a;
+  Vector_3 bd=d-b;
+  Vector_3 cd=d-c;
+
+  assert(orientation(ad, bd, cd) == POSITIVE);
+
+  assert(orientation(a,b,c) != POSITIVE); // The origin is below abc
+  bool ori_adb=orientation(a,d,b) == POSITIVE;
+  bool ori_bdc=orientation(b,d,c) == POSITIVE;
+  bool ori_cda=orientation(c,d,a) == POSITIVE;
+
+  //If the origin is inside the tetrahedron, no direction to return
+  if(!ori_adb && !ori_bdc && !ori_cda)
+    return NULL_VECTOR;
+
+  //Compute the position of the origin with adb
+  Vector_3 n_adb =cross(ad, bd);
+  bool on_left_of_ad = strictly_on_left_of_triangle_edge(-d, n_adb, ad, k);
+  bool on_right_of_bd = strictly_on_left_of_triangle_edge(-d, n_adb, -bd, k);
+
+  if(ori_adb && !on_left_of_ad && !on_right_of_bd){
+    // The origin is above the triangle adb
+    simplex[2]=simplex[3];
+    simplex.pop_back();
+    return n_adb;
+  }
+
+  //Compute the position of the origin with cda
+  Vector_3 n_cda =cross(cd, ad);
+  bool on_left_of_cd = strictly_on_left_of_triangle_edge(-d, n_cda, cd, k);
+  bool on_right_of_ad = strictly_on_left_of_triangle_edge(-d, n_cda, -ad, k);
+
+  if(ori_cda && !on_left_of_cd && !on_right_of_ad){
+    // The origin is above the triangle cda
+    simplex[1]=simplex[3];
+    simplex.pop_back();
+    return n_cda;
+  }
+
+  //Compute the position of the origin with bdc
+  Vector_3 n_bdc =cross(bd, cd);
+  bool on_left_of_bd = strictly_on_left_of_triangle_edge(-d, n_bdc, bd, k);
+  bool on_right_of_cd = strictly_on_left_of_triangle_edge(-d, n_bdc, -cd, k);
+
+  if(ori_bdc && !on_left_of_bd && !on_right_of_cd){
+    // The origin is above the triangle bdc
+    simplex[0]=simplex[3];
+    simplex.pop_back();
+    return n_bdc;
+  }
+
+  // The origin is not above a face
+  // Using the predicates already computed, we identify the edge below the origin
+  if(on_left_of_ad && on_right_of_ad){
+    if(is_negative(dot(ad, d))){
+      //The origin is above the vertex d
+      simplex[0]=simplex[3];
+      simplex.pop_back();
+      simplex.pop_back();
+      simplex.pop_back();
+      return -simplex[0];
+    }
+    Vector_3 dir=-cross(cross(ad, d), ad);
+    simplex[1]=simplex[3];
+    simplex.pop_back();
+    simplex.pop_back();
+    return dir;
+  }
+
+  if(on_left_of_bd && on_right_of_bd){
+    if(is_negative(dot(bd, d))){
+      //The origin is above the vertex d
+      simplex[0]=simplex[3];
+      simplex.pop_back();
+      simplex.pop_back();
+      simplex.pop_back();
+      return -simplex[0];
+    }
+    Vector_3 dir=-cross(cross(bd, d), bd);
+    simplex[0]=simplex[3];
+    simplex.pop_back();
+    simplex.pop_back();
+    return dir;
+  }
+
+  assert(on_left_of_cd && on_right_of_cd);
+  if(is_negative(dot(cd, d))){
+    //The origin is above the vertex d
+    simplex[0]=simplex[3];
+    simplex.pop_back();
+    simplex.pop_back();
+    simplex.pop_back();
+    return -simplex[0];
+  }
+  Vector_3 dir=-cross(cross(cd, d), cd);
+  simplex[0]=simplex[2];
+  simplex[1]=simplex[3];
+  simplex.pop_back();
+  simplex.pop_back();
+  return dir;
+}
+
+template<typename K, typename Vector_3>
+Vector_3 dir_to_origin(boost::container::small_vector<Vector_3, 4>& simplex, const K &k) {
+  const typename K::Construct_cross_product_vector_3 cross=k.construct_cross_product_vector_3_object();
+  const typename K::Compute_scalar_product_3 dot=k.compute_scalar_product_3_object();
+
+  switch(simplex.size()){
+  case 0:
+    return Vector_3(1,0,0);
+    break;
+  case 1:
+    return -simplex[0];
+    break;
+  case 2:
+  {
+    Vector_3 &a = simplex[0];
+    Vector_3 &b = simplex[1];
+    Vector_3 ab=b-a;
+    typename K::FT ratio=dot(ab, -a);
+    if(is_negative(ratio)){
+      simplex.pop_back();
+      return -a;
+    }
+    if(ratio>ab.squared_length()){
+      Vector_3 dir=-b;
+      simplex[0]=simplex[1];
+      simplex.pop_back();
+      return dir;
+    }
+    return -cross(cross(ab, a), ab);
+    break;
+  }
+  case 3:
+    return triangle_dir_to_origin(simplex, k);
+    break;
+  default: //4
+    return tetrahedron_dir_to_origin(simplex, k);
+    break;
+  }
+}
+
+template<typename OK, typename K>
+struct Separation_distance_functor{
+
+  template<typename Convex, typename NamedParameters1, typename NamedParameters2>
+  typename OK::FT operator()(const Convex &a, const Convex &b, const NamedParameters1 &np1, const NamedParameters2 &np2) const{
+    using Vector_3= typename K::Vector_3;
+    using FT = typename K::FT;
+
+    const int INTER_MAX_ITER=0;
+    boost::container::small_vector<Vector_3, 4> simplex;
+
+    Cartesian_converter<K, OK> converter;
+
+    //The Kernel used is EPECK and not the one of the named parameter
+    const K &k=K();
+
+    unsigned long planeStatPerPair = 0;
+    do {
+#ifdef CGAL_CONVEX_HULL_3_DISTANCE_VERBOSE
+      std::cout << "\nIteration " << planeStatPerPair << std::endl;
+      std::cout << "Simplex size: " << simplex.size() << ", points of the simplex:"
+      for(const Vector_3 &v: simplex)
+        std::cout << v << std::endl;
+#endif
+
+      Vector_3 dir=dir_to_origin<K>(simplex, k);
+      if(dir==NULL_VECTOR) return 0;
+
+      Vector_3 sp = Convex_hull_3::internal::extreme_point_3_wrapper(a, dir.direction(), np1) - Convex_hull_3::internal::extreme_point_3_wrapper(b, -dir.direction(), np2);
+      if(sp==NULL_VECTOR) return 0;
+      if(INTER_MAX_ITER!=0 && (++planeStatPerPair >= INTER_MAX_ITER)) return 0;
+
+#ifdef CGAL_CONVEX_HULL_3_DISTANCE_VERBOSE
+      std::cout << "Direction to origin: " << dir << std::endl;
+      std::cout << "Support point: " << sp << std::endl;
+#endif
+
+      simplex.push_back(sp);
+      //If the simplex is degenerated, the algorithm has terminate
+      if( (simplex.size()==4 && orientation(ORIGIN+simplex[0],ORIGIN+simplex[1],ORIGIN+simplex[2],ORIGIN+simplex[3])==COPLANAR) ||
+          (simplex.size()==3 && collinear(ORIGIN+simplex[0], ORIGIN+simplex[1], ORIGIN+simplex[2])) ||
+          (simplex.size()==2 && (simplex[0]==simplex[1])) )
+      {
+        simplex.pop_back();
+        break;
+      }
+
+      assert(++planeStatPerPair<=30);
+    } while( true );
+    FT dist;
+    typename K::Point_3 ORIGIN(0,0,0);
+    if(simplex.size()>=3)
+      dist = squared_distance(typename K::Plane_3(ORIGIN+simplex[0], ORIGIN+simplex[1], ORIGIN+simplex[2]), ORIGIN);
+    else if(simplex.size()==2)
+      dist = squared_distance(typename K::Line_3(ORIGIN+simplex[0], ORIGIN+simplex[1]), ORIGIN);
+    else
+      dist = simplex[0].squared_length();
+    return converter(dist);
+  }
+};
+
+} // end of predicates_impl namespace
+
+
+/**
+* \ingroup PkgConvexHull3Intersections
+*
+* \brief compute the separation distance of two convex hulls returning zero if they intersect.
+*
+* Input hulls can be provided as a range of points or as a graph, and may be of different types.
+* Furthermore, when many intersection queries use the same object, one should wrap that input in the class `CGAL::Convex_hull_hierarchy_3` as to construct an optimized view of the convex hull and accelerate intersection tests.
+* This is especially true when the convex hull is made of a large number of vertices (see `CGAL::Convex_hull_hierarchy_3` for more details).
+*
+* @tparam Convex_1 is one of the following types:\n
+*  - a model of `ConstRange`
+*  - a model of `VertexListGraph` and `AdjacencyGraph`
+*  - an instance of `CGAL::Convex_hull_hierarchy_3`
+* @tparam Convex_2 same as `Convex_1`
+* @tparam NamedParameters_1 a sequence of \ref bgl_namedparameters "Named Parameters"
+* @tparam NamedParameters_2 a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param ch1 the first convex hull
+* @param ch2 the second convex hull
+* @param np1 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* @param np2 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{when `ch1` (`ch2`) is a range, it is a property map associating points to its elements}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose value types are the same for `ch1` and `ch2`}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*     \cgalParamExtra{used only if `ch1` (`ch2`) is model of `ConstRange`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{when `ch1` (`ch2`) is a mesh, it is a property map associating points to its vertices}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose value types are the same for `ch1` and `ch2`}
+*     \cgalParamDefault{boost::get(CGAL::vertex_point, g)}
+*     \cgalParamExtra{used only if `ch1` (`ch2`) is model of `VertexListGraph` and `AdjacencyGraph`.}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` must be available in `Convex_1` (`Convex_2`).}
+*   \cgalParamNEnd
+*
+*   \cond SKIP_IN_MANUAL
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{An instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal kernel deduced from the point type of the input, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{`np1` only}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{number_of_iterations}
+*     \cgalParamDescription{if not `0` (no limit), indicates the maximum number of iterations performed by the algorithm.
+*                           If this value is not `0`, then an intersection might be reported even if the convex hulls does not intersect.
+*                           However, if the convex hulls are reported not to intersect, this is guaranteed.}
+*     \cgalParamType{a positive integer convertible to `std::size_t`}
+*     \cgalParamExtra{`np1` only}
+*     \cgalParamDefault{`0`}
+*   \cgalParamNEnd
+*
+* \endcond
+*
+* \cgalNamedParamsEnd
+*
+* \see `CGAL::Convex_hull_hierarchy_3`
+*/
+template <class Convex_1, class Convex_2,
+          class NamedParameters_1 = parameters::Default_named_parameters,
+          class NamedParameters_2 = parameters::Default_named_parameters>
+typename internal::GetGeomTraitsFromConvex<Convex_1, NamedParameters_1>::type::FT
+separation_distance(const Convex_1& c1, const Convex_2& c2,
+                    const NamedParameters_1& np1 = parameters::default_values(),
+                    const NamedParameters_2& np2 = parameters::default_values()){
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  //The function need exact computation to works correctly
+  using EPECK=Exact_predicates_exact_constructions_kernel;
+  using GT= typename internal::GetGeomTraitsFromConvex<Convex_1, NamedParameters_1>::type;
+  // GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+  return predicates_impl::Separation_distance_functor<GT, EPECK>()(c1, c2, np1, np2);
+}
+
+} // end of experimental
+
+}} // CGAL::Convex_hull_3 namespace
+
+#endif // CGAL_CONVEX_HULL_3_DISTANCE_H
diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/do_intersect.h b/Convex_hull_3/include/CGAL/Convex_hull_3/do_intersect.h
new file mode 100644
index 000000000000..824aaca5cd1d
--- /dev/null
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/do_intersect.h
@@ -0,0 +1,346 @@
+// Copyright (c) 2022 INRIA (France) and 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)     : Samuel Hornus, Sébastien Loriot and Léo Valque
+//
+
+#ifndef CGAL_CONVEX_HULL_3_DO_INTERSECT_H
+#define CGAL_CONVEX_HULL_3_DO_INTERSECT_H
+
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/Kernel_23/internal/Has_boolean_tags.h>
+
+#include <CGAL/Convex_hull_hierarchy_3.h>
+#include <CGAL/extreme_point_3.h>
+#include <CGAL/intersections.h>
+
+#include <CGAL/Convex_hull_3/internal/helpers.h>
+
+#include <CGAL/Container_helper.h>
+#include <CGAL/Dynamic_property_map.h>
+#include <CGAL/property_map.h>
+
+#include <boost/range/value_type.hpp>
+#include <boost/range/reference.hpp>
+
+#include <boost/graph/adjacency_list.hpp>
+
+#include <vector>
+
+namespace CGAL {
+
+namespace Convex_hull_3 {
+
+namespace internal {
+
+template <class Vector_3>
+Vector_3 LInf_normalize(const Vector_3 &v){
+  if(v==NULL_VECTOR)
+    return v;
+  auto l=(max)(abs(v.x()), (max)(abs(v.y()),abs(v.z())));
+  return v/l;
+}
+
+template <class Vector_3>
+struct SphericalPolygonElement {
+  Vector_3 vertex_; // A vertex of the spherical polygon
+  Vector_3 north_; // The north pole of the equatorial arc/edge leading OUT OF that vertex_ (arcs are oriented west-to-east, or CCW in a right-handed frame.
+  // In the spherical polygon (v0, n0), (v1, n1), (v2, n2), ... we have
+  // v1 = cross(n0, n1),  more generally: v_{i+1} = cross(n_i, n_{i+1})  and
+  // n1 = cross(v1, v2),  more generally: n_i     = cross(v_i, v_{i+1}).
+  SphericalPolygonElement(){}
+  SphericalPolygonElement(const Vector_3 & n) : north_(n) {}
+  SphericalPolygonElement(const Vector_3 & v, const Vector_3 & n) : vertex_(v), north_(n) {}
+};
+
+template <class Vector_3>
+struct SphericalPolygon : public std::vector<SphericalPolygonElement<Vector_3>> {
+
+  typedef std::vector<SphericalPolygonElement<Vector_3>> Base;
+  typedef typename Base::iterator iterator;
+  typedef typename Base::const_iterator const_iterator;
+  typedef typename Kernel_traits<Vector_3>::Kernel Kernel;
+  typedef typename Kernel::FT NT;
+
+  SphericalPolygon() {
+    this->reserve(16);
+  }
+
+  Vector_3 averageDirection() const {
+    // PRECONDITION : all vertices (that are unit vectors on the sphere) are normalized.
+    switch( this->size() ) {
+      case 0 : return Vector_3(1,0,0); break; // An arbitrary one
+      case 1 : return this->begin()->north_; break;
+      case 2 :{   // The two vertices are opposite so we do not take their mean
+                  // We want a direction with negative dot with bith north_
+                  // We take the two perpendicular vector of resp north_0, and north_1 that lies on the same circle on them
+                  // And we take a barycenter of them
+                  Vector_3 perp1= LInf_normalize(cross_product((*this)[0].north_, (*this)[0].vertex_));
+                  Vector_3 perp2= LInf_normalize(cross_product((*this)[1].north_, (*this)[1].vertex_));
+                  CGAL_assertion(is_positive((perp1+perp2) * (*this)[0].north_));
+                  CGAL_assertion(is_positive((perp1+perp2) * (*this)[1].north_));
+                  return (perp1 + perp2);
+                  // return (*this)[0].north_ + (*this)[1].north_; //Old, need that both north_ was L2 normalized
+               } break;
+      default : {
+                  Vector_3 avg(NULL_VECTOR);
+                  for( const SphericalPolygonElement<Vector_3> & v : *this )
+                    avg += v.vertex_;
+                  return avg;
+                } break;
+    }
+  }
+
+  void clip(const Vector_3& clipNorth, SphericalPolygon<Vector_3> &result) const {
+    const int n = this->size();
+    result.clear();
+    switch( n ) {
+      case 0 : {
+                  result.emplace_back(clipNorth);
+                  break;
+               }
+      case 1 : {
+                 result = (*this);
+                //  NT dot = this->begin()->north_ * clipNorth;
+                 Vector_3 v(LInf_normalize(cross_product(clipNorth, this->begin()->north_)));
+                 if(v==NULL_VECTOR /* && is_negative(dot) */){
+                  CGAL_assertion(is_negative(clipNorth * this->begin()->north_)); //Could not be two equal hemispheres
+                   // intersection of two opposite hemispheres ==> empty
+                   result.clear();
+                   break;
+                 }
+
+                 result.begin()->vertex_ = v;
+                 result.emplace_back(-v, clipNorth);
+                 break;
+               }
+      case 2 : {
+                 result = (*this);
+                 iterator next = result.begin();
+                 iterator cur = next++;
+                 NT vDot = this->begin()->vertex_ * clipNorth;
+                 if( is_positive(vDot) ) {
+                   // we'll get a triangle
+                   next->vertex_ = LInf_normalize(cross_product(clipNorth, next->north_));
+                   Vector_3 v( LInf_normalize(cross_product(cur->north_, clipNorth)));
+                   result.emplace(next, v, clipNorth);
+                 } else if( is_negative(vDot) ) {
+                   // we'll get a triangle
+                   cur->vertex_ =  LInf_normalize(cross_product(clipNorth, cur->north_));
+                   Vector_3 v( LInf_normalize(cross_product(next->north_, clipNorth)));
+                   result.emplace_back(v, clipNorth);
+                 } else {
+                   // we keep a moon crescent
+                   NT curTest(clipNorth *  cross_product(cur->north_, cur->vertex_));
+                   Vector_3 nextTest( cross_product(next->north_, next->vertex_));
+                   if( is_positive(curTest) ) {
+                    //The clipNorth is not between the two previous ones
+                    CGAL_assertion(!is_positive(clipNorth * nextTest));
+                     next->north_ = clipNorth;
+                     cur->vertex_ =  LInf_normalize(cross_product(next->north_, cur->north_));
+                     next->vertex_ = -cur->vertex_;
+                   } else {
+                     if( is_positive(clipNorth * nextTest) ) {
+                       cur->north_ = clipNorth;
+                       next->vertex_ =  LInf_normalize(cross_product(cur->north_, next->north_));
+                       cur->vertex_ = -next->vertex_;
+                     } else {
+                       //killed the crescent
+                       result.clear();
+                     }
+                   }
+                 }
+                 break;
+               }
+      default : { // n >= 3
+                  int nbKept(0);
+                  const_iterator cur = this->begin();
+                  NT nextDot, curDot = clipNorth * cur->vertex_;
+                  while( cur != this->end() ) {
+                    if( cur+1 == this->end() )
+                      nextDot = clipNorth * this->begin()->vertex_;
+                    else
+                      nextDot = clipNorth * (cur+1)->vertex_;
+                    if( is_positive(curDot) ) { // cur is "IN"
+                      ++nbKept;
+                      result.push_back(*cur);
+                      if( is_negative(nextDot) ) { // next is "OUT"
+                        result.emplace_back( LInf_normalize(cross_product(cur->north_, clipNorth)), clipNorth);
+                      }
+                    } else if( !is_negative(curDot) ) { // cur is "ON" the clipping plane
+                      ++nbKept;
+                      if ( is_negative(nextDot) ) // next is "OUT"
+                        result.emplace_back(cur->vertex_, clipNorth);
+                      else
+                        result.push_back(*cur);
+                    } else { // cur is "OUT"
+                      if ( is_positive(nextDot) ) { // next is "IN"
+                        result.emplace_back( LInf_normalize(cross_product(clipNorth, cur->north_)), cur->north_);
+                      }
+                    }
+                    curDot = nextDot;
+                    ++cur;
+                  }
+                  if( (result.size() < 3/*too small*/) || ((nbKept == n)/*no change*/) ) {
+                    result.clear();
+                  }
+                  break;
+                }
+    }
+  }
+};
+
+template<typename K>
+struct Functor_do_intersect{
+
+  template<typename Convex, typename NamedParameter1, typename NamedParameter2>
+  bool operator()(const Convex &a, const Convex &b, const NamedParameter1 &np1, const NamedParameter2 &np2) const{
+    using Vector_3= typename K::Vector_3;
+    const int INTER_MAX_ITER=0;
+    SphericalPolygon<Vector_3> positiveBound, tempPoly;
+    positiveBound.clear();
+    unsigned long planeStatPerPair = 0;
+    do {
+      Vector_3 dir = positiveBound.averageDirection();
+      Vector_3 sp = extreme_point_3_wrapper(a, dir.direction(), np1) - extreme_point_3_wrapper(b, -dir.direction(), np2);
+      if(sp==NULL_VECTOR) return true;
+      if(is_negative(sp * dir)) return false;
+      if(INTER_MAX_ITER!=0 && (++planeStatPerPair >= INTER_MAX_ITER)) return true;
+      positiveBound.clip(-sp, tempPoly); positiveBound.swap(tempPoly);
+    } while( !positiveBound.empty() );
+    return true;
+  }
+};
+
+} // end of internal namespace
+
+//Do_intersect_traits, only used to filter the kernel
+template<typename K,
+         typename IK=K,
+         typename Converter=Cartesian_converter<K, K>,
+         bool Has_filtered_predicates_ = CGAL::internal::Has_filtered_predicates<K>::value>
+struct Do_intersect_traits;
+
+template<typename K, typename IK, typename Converter>
+struct Do_intersect_traits<K, IK, Converter, false>{
+  typedef internal::Functor_do_intersect<K> Do_intersect;
+  Do_intersect do_intersect_object() const {
+    return Do_intersect();
+  }
+};
+template<typename K, typename Converter>
+struct Do_intersect_traits<K, K, Converter, true> {
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::Exact_kernel::Vector_3 EVector_3;
+  typedef typename K::Approximate_kernel::Vector_3 FVector_3;
+
+  typedef Cartesian_converter<K, K>  IdentityConverter;
+  typedef typename K::C2E C2E;
+  typedef typename K::C2F C2F;
+
+  typedef Do_intersect_traits<typename K::Exact_kernel, K, C2E> Exact_traits;
+  typedef Do_intersect_traits<typename K::Approximate_kernel, K, C2F> Filtering_traits;
+
+  //The conversion are made lazely by Do_intersect, we thus use IdentityConverter in Filtered_predicate
+  typedef Filtered_predicate<
+              typename Exact_traits::Do_intersect,
+              typename Filtering_traits::Do_intersect,
+              IdentityConverter,
+              IdentityConverter>  Do_intersect;
+
+  Do_intersect do_intersect_object() const
+  {
+    typename Exact_traits::Do_intersect pe = Exact_traits().do_intersect_object();
+    typename Filtering_traits::Do_intersect pf = Filtering_traits().do_intersect_object();
+
+    return Do_intersect(pe, pf);
+  }
+};
+
+/**
+* \ingroup PkgConvexHull3Intersections
+*
+* \brief checks whether two convex hulls intersect or not.
+*
+* Input hulls can be provided as a range of points or as a graph, and may be of different types.
+* Furthermore, when many intersection queries use the same object, one should wrap that input in the class `CGAL::Convex_hull_hierarchy_3` as to construct an optimized view of the convex hull and accelerate intersection tests.
+* This is especially true when the convex hull is made of a large number of vertices (see `CGAL::Convex_hull_hierarchy_3` for more details).
+*
+* @tparam Convex_1 is one of the following types:\n
+*  - a model of `ConstRange`
+*  - a model of `VertexListGraph` and `AdjacencyGraph`
+*  - an instance of `CGAL::Convex_hull_hierarchy_3`
+* @tparam Convex_2 same as `Convex_1`
+* @tparam NamedParameters_1 a sequence of \ref bgl_namedparameters "Named Parameters"
+* @tparam NamedParameters_2 a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param ch1 the first convex hull
+* @param ch2 the second convex hull
+* @param np1 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* @param np2 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{when `ch1` (`ch2`) is a range, it is a property map associating points to its elements}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose value types are the same for `ch1` and `ch2`}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*     \cgalParamExtra{used only if `ch1` (`ch2`) is model of `ConstRange`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{when `ch1` (`ch2`) is a mesh, it is a property map associating points to its vertices}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose value types are the same for `ch1` and `ch2`}
+*     \cgalParamDefault{boost::get(CGAL::vertex_point, g)}
+*     \cgalParamExtra{used only if `ch1` (`ch2`) is model of `VertexListGraph` and `AdjacencyGraph`.}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` must be available in `Convex_1` (`Convex_2`).}
+*   \cgalParamNEnd
+*
+*   \cond SKIP_IN_MANUAL
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{An instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal kernel deduced from the point type of the input, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{`np1` only}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{number_of_iterations}
+*     \cgalParamDescription{if not `0` (no limit), indicates the maximum number of iterations performed by the algorithm.
+*                           If this value is not `0`, then an intersection might be reported even if the convex hulls does not intersect.
+*                           However, if the convex hulls are reported not to intersect, this is guaranteed.}
+*     \cgalParamType{a positive integer convertible to `std::size_t`}
+*     \cgalParamExtra{`np1` only}
+*     \cgalParamDefault{`0`}
+*   \cgalParamNEnd
+*
+* \endcond
+*
+* \cgalNamedParamsEnd
+*
+* \see `CGAL::Convex_hull_hierarchy_3`
+*/
+template <class Convex_1, class Convex_2,
+          class NamedParameters_1 = parameters::Default_named_parameters,
+          class NamedParameters_2 = parameters::Default_named_parameters>
+bool do_intersect(const Convex_1& ch1, const Convex_2& ch2,
+                  const NamedParameters_1& np1 = parameters::default_values(),
+                  const NamedParameters_2& np2 = parameters::default_values()){
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+  using GT= typename internal::GetGeomTraitsFromConvex<Convex_1, NamedParameters_1>::type;
+  // GT gt = choose_parameter<GT>(get_parameter(np1, internal_np::geom_traits));
+  return Do_intersect_traits<GT>().do_intersect_object()(ch1, ch2, np1, np2);
+}
+
+}} // CGAL::Convex_hull_3 namespace
+
+#endif // CGAL_CONVEX_HULL_3_DO_INTERSECT_H
diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/internal/helpers.h b/Convex_hull_3/include/CGAL/Convex_hull_3/internal/helpers.h
new file mode 100644
index 000000000000..4b926757c6e9
--- /dev/null
+++ b/Convex_hull_3/include/CGAL/Convex_hull_3/internal/helpers.h
@@ -0,0 +1,110 @@
+// Copyright (c) 2025 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_CONVEX_HULL_INTERNAL_HELPERS_H
+#define CGAL_CONVEX_HULL_INTERNAL_HELPERS_H
+
+#include <CGAL/IO/helpers.h>
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/Container_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/Convex_hull_hierarchy_3.h>
+#include <CGAL/extreme_point_3.h>
+
+#ifdef CGAL_PROFILE_CONVEX_HULL_DO_INTERSECT
+      size_t nb_visited=0;
+#endif
+
+namespace CGAL::Convex_hull_3::internal{
+
+template <class T>
+inline constexpr bool is_instance_of_CHH = false;
+
+template <class Arg>
+inline constexpr bool is_instance_of_CHH< Convex_hull_hierarchy_3<Arg> > = true;
+
+// template class to deduce the GT from Convex and NamedParameters
+template<class Convex,
+         class NamedParameters,
+         bool Is_range=CGAL::IO::internal::is_Range_v<Convex> >
+struct GetGeomTraitsFromConvex{
+  typedef typename GetGeomTraits<Convex, NamedParameters>::type type;
+};
+
+template<class Convex, class NamedParameters>
+struct GetGeomTraitsFromConvex<Convex, NamedParameters, true>{
+  typedef typename Point_set_processing_3_np_helper<Convex, NamedParameters>::Geom_traits type;
+};
+
+template<class Mesh, class NamedParameters>
+struct GetGeomTraitsFromConvex<Convex_hull_hierarchy_3<Mesh>, NamedParameters, false>{
+  typedef typename GetGeomTraits<Mesh, NamedParameters>::type type;
+};
+
+template <class Convex, class Direction_3, class NamedParameters>
+typename Kernel_traits<Direction_3>::Kernel::Point_3 extreme_point_3_wrapper(const Convex &c, const Direction_3 &dir, const NamedParameters &np){
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  if constexpr(Convex_hull_3::internal::is_instance_of_CHH<Convex>){
+    return c.extreme_point_3(dir);
+  } else if constexpr(CGAL::IO::internal::is_Range_v<Convex>){
+    using NP_helper= Point_set_processing_3_np_helper<Convex, NamedParameters>;
+
+    using Point_GT = typename NP_helper::Geom_traits;
+    using Default_GT = typename Kernel_traits<Direction_3>::Kernel;
+    using GT=typename internal_np::Lookup_named_param_def <
+        internal_np::geom_traits_t,
+        NamedParameters,
+        Default_GT
+      > ::type;
+
+    using Default_geom_traits_converter = Cartesian_converter<Point_GT, GT>;
+    using GTC=typename internal_np::Lookup_named_param_def <
+        internal_np::geom_traits_converter_t,
+        NamedParameters,
+        Default_geom_traits_converter
+      > ::type;
+    GTC converter = choose_parameter<GTC>(get_parameter(np, internal_np::geom_traits_converter));
+
+    auto point_map = NP_helper::get_const_point_map(c, np);
+    return converter(get(point_map, extreme_point_3(c, dir, np)));
+  } else { // Convex is a Graph
+    using GetGeomTraits = GetGeomTraits<Convex, NamedParameters>;
+    using GraphGT= typename GetGeomTraits::type;
+
+    using Default_GT = typename Kernel_traits<Direction_3>::Kernel;
+    using GT=typename internal_np::Lookup_named_param_def <
+      internal_np::geom_traits_t,
+      NamedParameters,
+      Default_GT
+    > ::type;
+
+    using Default_geom_traits_converter = Cartesian_converter<GraphGT, GT>;
+    using GTC=typename internal_np::Lookup_named_param_def <
+        internal_np::geom_traits_converter_t,
+        NamedParameters,
+        Default_geom_traits_converter
+      > ::type;
+    GTC converter = choose_parameter<GTC>(get_parameter(np, internal_np::geom_traits_converter));
+    auto point_map = GetVertexPointMap<Convex, NamedParameters>::get_const_map(np, c);
+    return converter(get(point_map, extreme_vertex_3(c, dir, np)));
+  }
+}
+
+}
+
+#endif // CGAL_CONVEX_HULL_INTERNAL_HELPERS_H
\ No newline at end of file
diff --git a/Convex_hull_3/include/CGAL/Convex_hull_hierarchy_3.h b/Convex_hull_3/include/CGAL/Convex_hull_hierarchy_3.h
new file mode 100644
index 000000000000..0497d0f27115
--- /dev/null
+++ b/Convex_hull_3/include/CGAL/Convex_hull_hierarchy_3.h
@@ -0,0 +1,341 @@
+// Copyright (c) 2025 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_CONVEX_HULL_HIERARCHY_3_H
+#define CGAL_CONVEX_HULL_HIERARCHY_3_H
+
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/property_map.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/Convex_hull_traits_3.h>
+#include <CGAL/convex_hull_3.h>
+#include <CGAL/extreme_point_3.h>
+#include <vector>
+
+namespace CGAL{
+
+  /// \ingroup PkgConvexHull3Intersections
+  /// \brief This class stores a convex hull in a data structure optimized for fast intersection tests.
+  ///
+  /// More specifically, the structure is optimized for `CGAL::extreme_vertex_3()`, which is used by `CGAL::Convex_hull_3::do_intersect()`.
+  /// The computational complexities of `CGAL::extreme_point_3()` and consequently `CGAL::Convex_hull_3::do_intersect()` is \f$O(n)\f$ for a range of points,
+  /// \f$O(\sqrt{n})\f$ on average for a graph, but only \f$O(\log{n})\f$ on average for a `Convex_hull_hierarchy_3`, where \f$n\f$ is the number of points of the object.
+  ///
+  /// Building this structure has linear complexity and is faster than computing a convex hull, but more costly than a single call to `CGAL::Convex_hull_3::do_intersect()`.
+  /// It is therefore relevant when many intersection queries are performed, particularly when a convex hull has a large number of vertices.
+  ///
+  /// To evaluate this class, we generated convex hulls by sampling \f$n\f$ random points on the unit sphere, computing their convex hulls, and performing intersection tests between them.
+  ///
+  /// Average performance measurements (times in milliseconds).
+  ///
+  /// | Sphere size | Hull build | Hierarchy build | do_intersect (Range) | do_intersect (Graph) | do_intersect (Hull Hierarchy) |
+  /// |-------------|-----------:|----------------:|---------------------:|---------------------:|------------------------------:|
+  /// | 40          | 0.06431  | 0.006908 | 0.00260 | 0.00174 | 0.00170  |
+  /// | 100         | 0.18069  | 0.016147 | 0.00594 | 0.00260 | 0.00272  |
+  /// | 350         | 0.68038  | 0.076661 | 0.01976 | 0.00735 | 0.00697  |
+  /// | 1,000       | 2.14852  | 0.207484 | 0.05664 | 0.02079 | 0.01292  |
+  /// | 3,500       | 7.49052  | 0.740683 | 0.19527 | 0.04876 | 0.01842  |
+  /// | 10,000      | 22.7057  | 2.22908  | 0.54977 | 0.07595 | 0.02789  |
+  /// | 35,000      | 93.7957  | 8.512    | 2.00917 | 0.13722 | 0.03730  |
+  /// | 100,000     | 399.377  | 32.0424  | 5.50728 | 0.39127 | 0.05841  |
+  ///
+  /// @tparam PolygonMesh The polygon mesh used to construct each level of the hierarchy. Must be a model of `VertexListGraph` and `MutableFaceGraph`.
+  ///         An internal property map for `CGAL::vertex_point_t` must be available, with a value type that is a `Point_3` typedef of a \cgal kernel.
+#if DOXYGEN_RUNNING
+template < class PolygonMesh>
+#else
+template < class PolygonMesh,
+           int MAX_HIERARCHY_DEPTH = 5,
+           int RATIO_BETWEEN_LEVELS = 24,
+           int MINSIZE_FOR_NEXT_LEVEL = RATIO_BETWEEN_LEVELS*6>
+#endif
+struct Convex_hull_hierarchy_3{
+  // parameterization of the hierarchy
+  /// @private
+  constexpr static size_t RATIO = RATIO_BETWEEN_LEVELS;
+  /// @private
+  constexpr static size_t MAXSIZE_FOR_NAIVE_SEARCH = RATIO;
+
+  /// @private
+  using VPM = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::const_type;
+
+  /// @private
+  using Point_3 = typename boost::property_traits<VPM>::value_type;
+
+  /// @private
+  using GT = typename Kernel_traits<Point_3>::Kernel;
+
+  /// @private
+  using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::vertex_descriptor;
+
+  /// @private
+  using V2VMap = typename boost::property_map<PolygonMesh, dynamic_vertex_property_t<vertex_descriptor> >::type;
+
+  /**
+  * Constructor taking a graph `g`. If `compute_convex_hull` is set to `true`, compute and use the convex hull of `g` instead.
+  *
+  * @tparam Graph a model of `VertexListGraph`
+  * @tparam NamedParameters a sequence of named parameters
+  *
+  * @param g the graph
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `g`}
+  *     \cgalParamType{a model of `ReadablePropertyMap` whose value type is a point type}
+  *     \cgalParamDefault{If this parameter is omitted, an internal property map for ` CGAL::vertex_point_t`  must be available in `VertexListGraph`}
+  *   \cgalParamNEnd
+  *   \cgalParamNBegin{compute_convex_hull}
+  *     \cgalParamDescription{If set to `true`, the convex hull of `g` is computed. Otherwise, `g` must be convex.}
+  *     \cgalParamType{Boolean}
+  *     \cgalParamDefault{true}
+  *   \cgalParamNEnd
+  *   \cgalParamNBegin{random_seed}
+  *     \cgalParamDescription{Define the seed used by the hierarchy.}
+  *     \cgalParamType{unsigned int}
+  *     \cgalParamDefault{The seed used by `CGAL::get_default_random()`.}
+  *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  */
+  template <typename Graph,
+            typename NamedParameters=parameters::Default_named_parameters>
+  Convex_hull_hierarchy_3(const Graph &g, const NamedParameters& np = parameters::default_values()){
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+    using parameters::is_default_parameter;
+
+    PolygonMesh ch;
+    bool compute_convex_hull = choose_parameter(get_parameter(np, internal_np::compute_convex_hull), true);
+    size_t seed = choose_parameter(get_parameter(np, internal_np::random_seed), 0);
+    Random rng = (is_default_parameter<NamedParameters, internal_np::random_seed_t>::value)? get_default_random(): Random(seed);
+    if(compute_convex_hull)
+      convex_hull_3(g, ch, np);
+    else
+      copy_face_graph(g, ch, np);
+    hierarchy_sm.reserve(MAX_HIERARCHY_DEPTH);
+    hierarchy_sm.push_back(std::move(ch));
+    init_hierarchy(rng);
+  };
+
+   /**
+  * Constructor taking the points in the range `[first, last)`.
+  *
+  * @tparam PointIterator is an input iterator with a value type equivalent to `Traits::Point_3`
+  * @tparam Traits is a model of the concept ConvexHullTraits_3. For the purposes of checking the postcondition that the convex hull is valid,
+  *         `Traits` must also be a model of the concept `IsStronglyConvexTraits_3`.
+  */
+  template <typename PointIterator, typename Traits = Convex_hull_traits_3<GT , PolygonMesh> >
+  Convex_hull_hierarchy_3(PointIterator begin, PointIterator end, const Traits &traits = Traits()){
+    PolygonMesh ch;
+    convex_hull_3(begin, end, ch, traits);
+    hierarchy_sm.reserve(MAX_HIERARCHY_DEPTH);
+    hierarchy_sm.push_back(std::move(ch));
+    init_hierarchy(get_default_random(), traits);
+  };
+
+  // /// Returns the index of the deepest level of the hierarchy (the one with the fewest vertices).
+  // /// This corresponds to the number of levels minus one.
+  /// @private
+  std::size_t maxlevel() const{
+    return hierarchy_sm.size()-1;
+  }
+
+  // /// returns the mesh of the corresponding level of the hierarchy.
+  /// @private
+  const PolygonMesh& mesh(std::size_t level) const{
+    return hierarchy_sm[level];
+  }
+
+  /// returns a const reference to the mesh stored by the class.
+  const PolygonMesh& mesh() const{
+    return hierarchy_sm[0];
+  }
+
+  /// @private
+  template <typename Direction>
+  typename Kernel_traits<Direction>::Kernel::Point_3 extreme_point_3(const Direction &dir) const {
+    using Dir_GT = typename Kernel_traits<Direction>::Kernel;
+    using GTC = Cartesian_converter<GT, Dir_GT>;
+    GTC converter;
+    VPM vpm = get_const_property_map(vertex_point, hierarchy_sm[0]);
+    return converter(get(vpm, extreme_vertex_3(dir)));
+  }
+
+  /**
+  * returns the vertex of the convex hull whose geometric position projected onto that direction is maximal.
+  * If not unique, a single vertex is returned.
+  *
+  * @tparam Direction must be a `Direction_3` typedef of a \cgal kernel
+  *
+  * @param dir the direction
+  */
+  template <typename Direction>
+  vertex_descriptor extreme_vertex_3(const Direction &dir) const
+  {
+    using Dir_GT = typename Kernel_traits<Direction>::Kernel;
+    using GTC = Cartesian_converter<GT, Dir_GT>;
+    Dir_GT gt;
+    GTC converter;
+
+    auto vector_3 = gt.construct_vector_3_object();
+    auto csp = gt.compare_scalar_product_3_object();
+
+    VPM vpm = get_const_property_map(vertex_point, hierarchy_sm[0]);
+    const typename Dir_GT::Vector_3 &vdir = dir.vector();
+
+    size_t level=maxlevel();
+
+    const PolygonMesh &sm = hierarchy_sm[level];
+    if(level == 0)
+      return CGAL::extreme_vertex_3(sm, dir);
+
+    vertex_descriptor argmax = *vertices(sm).begin();
+    Vector_3 vec_max = vector_3(ORIGIN, converter(get(vpm, get(to_base_maps[level], argmax))));
+    if(vertices(sm).size() <= MAXSIZE_FOR_NAIVE_SEARCH){
+      //If maxlevel is small, we simply go through all its vertices
+      for(auto vh=++(vertices(sm).begin()); vh!=vertices(sm).end(); ++vh){
+        vertex_descriptor v = *vh;
+        Vector_3 vec = vector_3(ORIGIN, converter(get(vpm, get(to_base_maps[level], v))));
+        if(csp(vec_max, vdir, vec, vdir)==SMALLER){
+          vec_max=vec;
+          argmax=v;
+        }
+      }
+
+      // Go to level below
+      if(level>0){
+        argmax= get(next_in_hierarchy_maps[level], argmax);
+        --level;
+      } else {
+        return argmax;
+      }
+    }
+
+    for(; true; --level){
+      // Starting from the vertex of the previous level, we walk on the graph
+      // along neighbors that increase the "score"
+      const PolygonMesh &csm = hierarchy_sm[level];
+      const V2VMap &cbase = to_base_maps[level];
+      const V2VMap &cnext = next_in_hierarchy_maps[level];
+      bool is_local_max;
+      do{
+        is_local_max=true;
+        for(vertex_descriptor v: vertices_around_target(argmax ,csm)){
+          Vector_3 vec = vector_3(ORIGIN, converter(get(vpm, get(cbase, v))));
+          if(csp(vec_max, vdir, vec, vdir)==SMALLER){
+            vec_max = vec;
+            argmax = v;
+            is_local_max = false; // repeat with the new vertex
+            break;
+          }
+        }
+      } while(!is_local_max);
+      if(level>0){
+        argmax= get(cnext, argmax);
+      } else {
+        return argmax;
+      }
+    }
+  }
+
+private:
+  template <typename Random, typename Traits = Convex_hull_traits_3<GT , PolygonMesh> >
+  void init_hierarchy(Random &rng = get_default_random(), const Traits &traits = Traits()){
+    VPM vpm = get_const_property_map(vertex_point, hierarchy_sm[0]);
+
+    size_t size=vertices(hierarchy_sm[0]).size();
+    size_t level=0;
+
+    to_base_maps.reserve(MAX_HIERARCHY_DEPTH);
+    next_in_hierarchy_maps.reserve(MAX_HIERARCHY_DEPTH);
+
+    V2VMap to_base_map = get(dynamic_vertex_property_t<vertex_descriptor>(), hierarchy_sm[0]);
+    for(vertex_descriptor v : vertices(hierarchy_sm[0]))
+      put(to_base_map, v, v);
+    next_in_hierarchy_maps.push_back(V2VMap());
+    to_base_maps.push_back(std::move(to_base_map));
+
+    while(size>MINSIZE_FOR_NEXT_LEVEL && level<MAX_HIERARCHY_DEPTH-1){
+      // Create a new level
+      std::vector<Point_3> select_points;
+      std::map<Point_3, vertex_descriptor> select_vertices;
+      select_points.reserve(2*size/RATIO);
+      PolygonMesh& above_sm=hierarchy_sm[level];
+      ++level;
+
+      // Select randomly vertices of the new level
+      while(select_points.size()==0)
+        for(vertex_descriptor v: vertices(above_sm))
+          if(rng.get_int(0,RATIO-1)==0){
+            vertex_descriptor base = get(to_base_maps[level-1], v);
+            select_points.emplace_back(get(vpm, base));
+            select_vertices[select_points.back()] = v;
+          }
+
+      // Compute the graph of the new level
+      PolygonMesh new_sm;
+      convex_hull_3(select_points.begin(), select_points.end(), new_sm, traits);
+      hierarchy_sm.push_back(std::move(new_sm));
+
+      auto temp_vpm = get_property_map(vertex_point, hierarchy_sm.back());
+
+      V2VMap to_base_map = get(CGAL::dynamic_vertex_property_t<vertex_descriptor>(),  hierarchy_sm.back());
+      V2VMap next_in_hierarchy_map = get(CGAL::dynamic_vertex_property_t<vertex_descriptor>(),  hierarchy_sm.back());
+
+      for(vertex_descriptor v : vertices( hierarchy_sm.back())){
+        vertex_descriptor next = select_vertices[get(temp_vpm, v)];
+        vertex_descriptor base = get(to_base_maps[level-1], next);
+        put(next_in_hierarchy_map, v, next);
+        put(to_base_map, v, base);
+      }
+
+      next_in_hierarchy_maps.push_back(std::move(next_in_hierarchy_map));
+      to_base_maps.push_back(std::move(to_base_map));
+      size=vertices( hierarchy_sm.back()).size();
+
+      CGAL_assertion(size!=0);
+    }
+  }
+
+  std::vector< PolygonMesh > hierarchy_sm;
+  // Given a vertex of level n, give the correspondant vertex of level n-1
+  std::vector<V2VMap> next_in_hierarchy_maps;
+  // Given a vertex of level n, give the correspondant vertex of level 0
+  std::vector<V2VMap> to_base_maps;
+};
+
+/**
+* \ingroup PkgConvexHull3Queries
+*
+* returns the vertex of the convex hull whose geometric position projected onto that direction is maximal.
+* If not unique, a single vertex is returned.
+*
+* @tparam Mesh a model of `VertexListGraph` and `MutableFaceGraph`
+* @tparam Direction is a model of `Kernel::Direction_3`.
+*
+* @param ch the convex hull
+* @param dir the direction
+*
+*/
+template <class PolygonMesh, class Direction>
+typename boost::graph_traits<PolygonMesh>::vertex_descriptor
+extreme_vertex_3(const Convex_hull_hierarchy_3<PolygonMesh> &ch, const Direction &dir)
+{
+  return ch.extreme_vertex_3(dir);
+}
+
+}
+
+#endif // CGAL_CONVEX_HULL_HIERARCHY_3_H
\ No newline at end of file
diff --git a/Convex_hull_3/include/CGAL/extreme_point_3.h b/Convex_hull_3/include/CGAL/extreme_point_3.h
new file mode 100644
index 000000000000..829a68375248
--- /dev/null
+++ b/Convex_hull_3/include/CGAL/extreme_point_3.h
@@ -0,0 +1,227 @@
+// Copyright (c) 2022 INRIA (France) and 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_EXTREME_POINT_3_H
+#define CGAL_EXTREME_POINT_3_H
+
+#include <CGAL/license/Convex_hull_3.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/Kernel_23/internal/Has_boolean_tags.h>
+
+#include <CGAL/Container_helper.h>
+#include <CGAL/Dynamic_property_map.h>
+#include <CGAL/property_map.h>
+
+#include <boost/range/value_type.hpp>
+#include <boost/range/reference.hpp>
+
+#include <boost/graph/adjacency_list.hpp>
+
+#include <vector>
+
+namespace CGAL {
+
+/**
+* \ingroup PkgConvexHull3Queries
+*
+* returns the point of the convex hull whose projection onto that direction is maximal.
+* If not unique, a single point is returned.
+*
+* @tparam PointRange is a model of `ConstRange`. The value type of its iterator is the key type of the named parameter `point_map`.
+* @tparam Direction must be a `Direction_3` typedef of a \cgal kernel
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param r the range of points
+* @param dir the direction
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of `r`}
+*     \cgalParamType{a model of `ReadablePropertyMap` with value type is a model `CGAL::Point_3`}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*   \cgalParamNEnd
+*
+* \cond SKIP_IN_MANUAL
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{An instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal kernel deduced from `Direction_3`, using `CGAL::Kernel_traits`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{geom_traits_converter}
+*     \cgalParamDescription{A Converter from the point type of `point_map` to the point type of `geom_traits`}
+*     \cgalParamType{a class model of `NT_Converter`}
+*     \cgalParamDefault{a \cgal `Cartesian_converter` deduced from `point_map` and `geom_traits`, using `CGAL::Kernel_traits`}
+*   \cgalParamNEnd
+*
+* \endcond
+*
+* \cgalNamedParamsEnd
+*
+* \return an instance of the key type of `point_map` parameter
+*/
+template <class PointRange, class Direction, class NamedParameters = parameters::Default_named_parameters>
+#if DOXYGEN_RUNNING
+Point_3
+#else
+typename Point_set_processing_3_np_helper<PointRange, NamedParameters>::Const_point_map::key_type
+#endif
+extreme_point_3(const PointRange& r, const Direction &dir, const NamedParameters &np = parameters::default_values()) {
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+  using NP_helper= Point_set_processing_3_np_helper<PointRange, NamedParameters>;
+
+  using PointMap= typename NP_helper::Const_point_map;
+  PointMap point_map = NP_helper::get_const_point_map(r, np);
+
+  using Point_GT = typename NP_helper::Geom_traits;
+  using Default_GT = typename Kernel_traits<Direction>::Kernel;
+  using GT=typename internal_np::Lookup_named_param_def <
+      internal_np::geom_traits_t,
+      NamedParameters,
+      Default_GT
+    > ::type;
+  GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+  using Vector_3 = typename GT::Vector_3;
+
+  using Default_geom_traits_converter = Cartesian_converter<Point_GT, GT>;
+  using GTC=typename internal_np::Lookup_named_param_def <
+      internal_np::geom_traits_converter_t,
+      NamedParameters,
+      Default_geom_traits_converter
+    > ::type;
+  GTC converter = choose_parameter<GTC>(get_parameter(np, internal_np::geom_traits_converter));
+
+  auto vector_3 = gt.construct_vector_3_object();
+  auto csp = gt.compare_scalar_product_3_object();
+
+  typename PointRange::const_iterator argmax=r.begin();
+  Vector_3 vec_max = vector_3(ORIGIN, converter(get(point_map, *argmax)));
+  for(typename PointRange::const_iterator it=++r.begin(); it!=r.end(); ++it){
+    Vector_3 vec = vector_3(ORIGIN, converter(get(point_map, *it)));
+    if(csp(vec_max, dir.vector(), vec, dir.vector())==SMALLER){
+      vec_max=vec;
+      argmax=it;
+    }
+  }
+  return *argmax;
+}
+
+/**
+* \ingroup PkgConvexHull3Queries
+*
+* returns the vertex of the convex hull whose geometric position projected onto that direction is maximal.
+* If not unique, a single vertex is returned.
+*
+* @tparam Graph is a model of `VertexListGraph` and `AdjacencyGraph`.
+* @tparam Direction must be a `Direction_3` typedef of a \cgal kernel
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param g the convex graph
+* @param dir the direction
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \pre The input graph must represent a convex object to guarantee a correct answer.
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{a property map associating points to the vertices of `g`}
+*     \cgalParamType{a model of `ReadablePropertyMap`}
+*     \cgalParamDefault{boost::get(CGAL::vertex_point, g)}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` must be available in `Graph`.}
+*   \cgalParamNEnd
+*
+* \cond SKIP_IN_MANUAL
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{An instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal kernel deduced from `Direction_3`, using `CGAL::Kernel_traits`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{geom_traits_converter}
+*     \cgalParamDescription{A converter from the point type of `vertex_point_map` to the point type of `geom_traits`}
+*     \cgalParamType{a class model of `NT_Converter`}
+*     \cgalParamDefault{a \cgal `Cartesian_converter` deduced from ` vertex_point_map`  and `geom_traits`, using `CGAL::Kernel_traits`}
+*   \cgalParamNEnd
+*
+* \endcond
+*
+* \cgalNamedParamsEnd
+*
+*/
+template <class Graph, class Direction, class NamedParameters = parameters::Default_named_parameters>
+typename boost::graph_traits<Graph>::vertex_descriptor
+extreme_vertex_3(const Graph& g, const Direction &dir, const NamedParameters &np = parameters::default_values()) {
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  using vertex_descriptor= typename boost::graph_traits<Graph>::vertex_descriptor;
+
+  using GetVertexPointMap = GetVertexPointMap<Graph, NamedParameters>;
+  using VPM = typename GetVertexPointMap::const_type;
+  VPM point_map = GetVertexPointMap::get_const_map(np, g);
+
+  using GetGeomTraits = GetGeomTraits<Graph, NamedParameters>;
+  using GraphGT= typename GetGeomTraits::type;
+
+  using Default_GT = typename Kernel_traits<Direction>::Kernel;
+  using GT=typename internal_np::Lookup_named_param_def <
+      internal_np::geom_traits_t,
+      NamedParameters,
+      Default_GT
+    > ::type;
+  GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+  using Vector_3 = typename GT::Vector_3;
+
+  using Default_geom_traits_converter = Cartesian_converter<GraphGT, GT>;
+  using GTC=typename internal_np::Lookup_named_param_def <
+      internal_np::geom_traits_converter_t,
+      NamedParameters,
+      Default_geom_traits_converter
+    > ::type;
+  GTC converter = choose_parameter<GTC>(get_parameter(np, internal_np::geom_traits_converter));
+
+  auto vector_3 = gt.construct_vector_3_object();
+  auto csp = gt.compare_scalar_product_3_object();
+
+  // If the number of vertices is small, simply test all vertices
+  if(vertices(g).size()<20)
+    return extreme_point_3(vertices(g), dir, parameters::point_map(point_map));
+
+  //Walks on the mesh to find a local maximun
+  vertex_descriptor argmax = *vertices(g).begin();
+  Vector_3 vec_max = vector_3(ORIGIN, converter(get(point_map,argmax)));
+  bool is_local_max;
+  do{
+    is_local_max=true;
+    for(auto v: vertices_around_target(argmax, g)){
+      Vector_3 vec = vector_3(ORIGIN, converter(get(point_map, v)));
+      if(csp(vec_max, dir.vector(), vec, dir.vector())==SMALLER){
+        vec_max = vec;
+        argmax = v;
+        is_local_max=false; // repeat with the new vertex
+        break;
+      }
+    }
+  }while(!is_local_max);
+  // Since convex, local maximum is a global maximum
+  return argmax;
+}
+
+} // CGAL namespace
+
+#endif // CGAL_EXTREME_POINT_3_H
diff --git a/Convex_hull_3/package_info/Convex_hull_3/dependencies b/Convex_hull_3/package_info/Convex_hull_3/dependencies
index bdd3aa15916c..4d6f9077648a 100644
--- a/Convex_hull_3/package_info/Convex_hull_3/dependencies
+++ b/Convex_hull_3/package_info/Convex_hull_3/dependencies
@@ -11,11 +11,13 @@ Distance_3
 Filtered_kernel
 HalfedgeDS
 Hash_map
+Homogeneous_kernel
 Installation
 Intersections_2
 Intersections_3
 Interval_support
 Kernel_23
+Kernel_d
 Modular_arithmetic
 Number_types
 Profiling_tools
diff --git a/Convex_hull_3/test/Convex_hull_3/ch3_test_do_intersect.cpp b/Convex_hull_3/test/Convex_hull_3/ch3_test_do_intersect.cpp
new file mode 100644
index 000000000000..239d2a8db9bb
--- /dev/null
+++ b/Convex_hull_3/test/Convex_hull_3/ch3_test_do_intersect.cpp
@@ -0,0 +1,309 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Simple_cartesian.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/Extreme_points_traits_adapter_3.h>
+
+#include <CGAL/Convex_hull_3/do_intersect.h>
+#include <CGAL/Convex_hull_3/distance.h>
+#include <CGAL/convex_hull_3.h>
+#include <CGAL/Convex_hull_hierarchy_3.h>
+
+#include <boost/property_map/vector_property_map.hpp>
+#include <CGAL/point_generators_3.h>
+
+#include <vector>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel        Epeck;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel      Epick;
+// typedef CGAL::Simple_cartesian<double>                           Cartesian;
+
+template<typename K>
+struct Sphere{
+  typedef typename K::FT FT;
+  typedef typename K::Point_3 P;
+
+  // do_intersect deduce Kernel from iterator type of the object, this an astuce for kernel deduction
+  typedef typename std::vector<P>::iterator iterator;
+
+  FT r;
+  P c;
+  Sphere(P c_, FT r_):r(r_), c(c_){}
+
+  template<class Vector_3, class Converter>
+  P extreme_point(const Vector_3 dir, const Converter& converter) const{
+    return converter(c)+converter(r)*dir/CGAL::sqrt(dir.squared_length());
+  }
+};
+template<typename K, class Vector_3, class Converter>
+typename K::Point_3 extreme_point(const Sphere<K> &sp, const Vector_3 &dir, const Converter& converter){
+  return sp.extreme_point(dir, converter);
+}
+
+template<typename K, typename Mesh>
+struct Test{
+  using P = typename K::Point_3;
+  using V = typename K::Vector_3;
+  using FT = typename K::FT;
+
+  typedef boost::vector_property_map<P> PMap;
+
+  void test(std::vector<P> &vec_a, std::vector<P> &vec_b, bool result, FT true_distance=FT(0)){
+    assert(CGAL::Convex_hull_3::do_intersect(vec_a, vec_b)==result);
+    assert(CGAL::Convex_hull_3::do_intersect(vec_b, vec_a)==result);
+
+    if(CGAL::is_zero(true_distance)){
+      assert((CGAL::Convex_hull_3::experimental::separation_distance(vec_a, vec_b)==0)==result);
+      assert((CGAL::Convex_hull_3::experimental::separation_distance(vec_b, vec_a)==0)==result);
+    } else {
+      assert(CGAL::Convex_hull_3::experimental::separation_distance(vec_a, vec_b)==true_distance);
+      assert(CGAL::Convex_hull_3::experimental::separation_distance(vec_b, vec_a)==true_distance);
+    }
+
+
+    Mesh sm_a, sm_b;
+    CGAL::convex_hull_3(vec_a.begin(), vec_a.end(), sm_a);
+    CGAL::convex_hull_3(vec_b.begin(), vec_b.end(), sm_b);
+    assert(CGAL::Convex_hull_3::do_intersect(sm_a, sm_b)==result);
+    assert(CGAL::Convex_hull_3::do_intersect(sm_b, sm_a)==result);
+
+    if(!CGAL::is_zero(true_distance)){
+      assert(CGAL::Convex_hull_3::experimental::separation_distance(sm_a, sm_b)==true_distance);
+      assert(CGAL::Convex_hull_3::experimental::separation_distance(sm_b, sm_a)==true_distance);
+    }
+
+    CGAL::Convex_hull_hierarchy_3<Mesh> hsm_a(vec_a.begin(),vec_a.end()), hsm_b(vec_b.begin(),vec_b.end());
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_a, hsm_b)==result);
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_b, hsm_a)==result);
+
+    CGAL::Convex_hull_hierarchy_3<Mesh> hsm_a_2(sm_a), hsm_b_2(sm_b);
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_a_2, hsm_b_2)==result);
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_b_2, hsm_a_2)==result);
+
+    // Test with Point map
+    PMap v_a, v_b; //Vectors de Point_3
+    for(size_t i=0; i<vec_a.size(); ++i)
+      v_a[i]=vec_a[i];
+    for(size_t i=0; i<vec_b.size(); ++i)
+      v_b[i]=vec_b[i];
+    std::vector< size_t > pm_a(vec_a.size(),0), pm_b(vec_b.size(), 0); //Vector of indexes
+    std::iota(pm_a.begin(), pm_a.end(), 0);
+    std::iota(pm_b.begin(), pm_b.end(), 0);
+    assert(CGAL::Convex_hull_3::do_intersect(pm_a, pm_b, CGAL::parameters::point_map(v_a), CGAL::parameters::point_map(v_b))==result);
+    assert(CGAL::Convex_hull_3::do_intersect(pm_b, pm_a, CGAL::parameters::point_map(v_b), CGAL::parameters::point_map(v_a))==result);
+
+    if constexpr(!std::is_same_v<Mesh, CGAL::Surface_mesh<P> >)
+      return;
+
+    using SM1 = CGAL::Surface_mesh<size_t>;
+
+    SM1 sm_a_pm, sm_b_pm;
+    CGAL::convex_hull_3(pm_a.begin(), pm_a.end(), sm_a_pm, CGAL::make_extreme_points_traits_adapter(v_a));
+    CGAL::convex_hull_3(pm_b.begin(), pm_b.end(), sm_b_pm, CGAL::make_extreme_points_traits_adapter(v_b));
+    auto vpm_a = make_compose_property_map(sm_a_pm.points(), v_a);
+    auto vpm_b = make_compose_property_map(sm_b_pm.points(), v_b);
+    assert(CGAL::Convex_hull_3::do_intersect(sm_a_pm, sm_b_pm, CGAL::parameters::vertex_point_map(vpm_a),
+                                                               CGAL::parameters::vertex_point_map(vpm_b))==result);
+    assert(CGAL::Convex_hull_3::do_intersect(sm_b_pm, sm_a_pm, CGAL::parameters::vertex_point_map(vpm_b),
+                                                               CGAL::parameters::vertex_point_map(vpm_a))==result);
+
+    CGAL::Convex_hull_hierarchy_3<Mesh> hsm_a_pm(sm_a_pm, CGAL::parameters::vertex_point_map(vpm_a));
+    CGAL::Convex_hull_hierarchy_3<Mesh> hsm_b_pm(sm_b_pm, CGAL::parameters::vertex_point_map(vpm_b));
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_a_pm, hsm_b_pm)==result);
+    assert(CGAL::Convex_hull_3::do_intersect(hsm_b_pm, hsm_a_pm)==result);
+  }
+
+  void test_cube()
+  {
+    std::vector<P> cube;
+    for(int x=0; x<2; ++x)
+      for(int y=0; y<2; ++y)
+        for(int z=0; z<2; ++z)
+            cube.push_back(P(x,y,z));
+
+    std::vector<P> inside(1, P(0.25,0.25,0.20));
+    std::vector<P> outside(1, P(-0.25,0.25,0.25));
+
+    test(cube, inside, true);
+    test(cube, outside, false, 1/16);
+
+    //Test intersection on vertex, edge, face
+    for(double x=0.; x<=1.; x+=0.5)
+      for(double y=0.; y<=1.; y+=0.5)
+        for(double z=0.; z<=1.; z+=0.5){
+          std::vector<P> vertex(1, P(x,y,z));
+          test(cube, vertex, true);
+        }
+
+    //Test between cubes
+    std::vector<P> cube_bis;
+    for(int x=0; x<2; ++x)
+      for(int y=0; y<2; ++y)
+        for(int z=0; z<2; ++z)
+            cube_bis.push_back(P(x,y,z));
+
+    auto transform=[](std::vector<P> &cube, const CGAL::Aff_transformation_3<K> &t){
+      for(auto &p: cube)
+        p=t(p);
+    };
+    //Test their intersection for many translations
+    for(double x=-1.5; x<=1.5; x+=0.5)
+      for(double y=-1.5; y<=1.5; y+=0.5)
+        for(double z=-1.5; z<=1.5; z+=0.5){
+          CGAL::Aff_transformation_3<K> t(CGAL::TRANSLATION, V(x,y,z));
+          transform(cube_bis, t);
+          test(cube, cube_bis, std::abs(x)<1.5 && std::abs(y)<1.5 && std::abs(z)<1.5);
+          transform(cube_bis, t.inverse());
+        }
+  }
+
+  void test_degenerate()
+  {
+    //Vertices
+    std::vector<P> origin(1, CGAL::ORIGIN);
+    std::vector<P> vertex1(1, P(1,0,0));
+    std::vector<P> vertex2(1, P(0,1,0));
+
+    test(origin, origin, true);
+    test(vertex1,vertex1, true);
+
+    test(origin, vertex1, false, 1);
+    test(vertex1, vertex2, false, 2);
+
+    //Segments
+    std::vector<P> seg1({P(0,0,0),P(2,0,0)});
+    std::vector<P> seg2({P(0,0,0),P(0,2,0)});
+    std::vector<P> seg3({P(0,2,0),P(2,2,0)});
+    std::vector<P> seg4({P(1,1,0),P(1,3,0)});
+
+    test(origin, seg1, true);
+    test(vertex1, seg1, true);
+
+    test(vertex2, seg1, false, 1);
+
+    test(seg1, seg1, true);
+    test(seg1, seg2, true);
+    test(seg3, seg4, true);
+
+    test(seg1, seg3, false, 4);
+    test(seg1, seg4, false, 1);
+
+    //Triangle
+    std::vector<P> tr1({P(0,0,0),P(2,0,0),P(0,2,0)});
+    std::vector<P> tr2({P(0,0,0),P(2,0,0),P(0,0,2)});
+    std::vector<P> tr3({P(0,0,0),P(2,0,2),P(0,0,2)});
+    std::vector<P> tr4({P(1,0,0),P(2,0,2),P(0,0,2)});
+    std::vector<P> tr5({P(0,0,2),P(2,0,2),P(0,2,2)});
+
+    test(tr1, tr1, true);
+    test(tr1, tr2, true);
+    test(tr1, tr3, true);
+    test(tr1, tr4, true);
+    test(tr3, tr4, true);
+
+    test(tr1, tr5, false, 4);
+    test(tr5, tr1, false, 4);
+  }
+
+  void test_half_sphere()
+  {
+    std::vector<P> half_sphere;
+    for(double phi=25./16.; phi>0; phi-=1./4.)
+      for(double theta=0; theta<2*CGAL_PI; theta+=0.25)
+        half_sphere.push_back(P(std::sin(phi) * std::cos(theta),
+                                std::sin(phi) * std::sin(theta),
+                                std::cos(phi)));
+
+    for(double x=-0.5; x<=0.5; x+=0.1)
+      for(double y=-0.5; y<=0.5; y+=0.1){
+        std::vector<P> outside(1, P(x,y,std::nextafter(std::cos(25./16),0)));
+        std::vector<P> inside(1, P(x,y,std::nextafter(std::cos(25./16),1)));
+
+        test(half_sphere, inside, true);
+        test(half_sphere, outside, false);
+      }
+  }
+
+  void test_random_tetrahedra(int N, CGAL::Random &r)
+  {
+    using Tet = typename K::Tetrahedron_3;
+    CGAL::Random_points_in_cube_3<P> gen (1, r);
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = *gen++;
+      P p1 = *gen++;
+      P p2 = *gen++;
+      P p3 = *gen++;
+
+      P q0 = *gen++;
+      P q1 = *gen++;
+      P q2 = *gen++;
+      P q3 = *gen++;
+
+      std::vector<P> a({p1,p2,p3,p0});
+      std::vector<P> b({q1,q2,q3,q0});
+      test(a, b, CGAL::do_intersect(Tet(p0, p1, p2, p3), Tet(q0, q1, q2, q3)));
+    }
+  }
+
+  void test_random_sphere(int N, int M, CGAL::Random &r)
+  {
+    auto random_sphere_point=[&](){
+      double a=r.get_double(0,2*CGAL_PI);
+      double b=r.get_double(-CGAL_PI/2,CGAL_PI/2);
+      return P( std::cos(a)*std::cos(b), std::sin(a)*std::cos(b), std::sin(b));
+    };
+    for(int i=0; i<M; ++i)
+    {
+      std::vector<P> sp_a;
+      std::vector<P> sp_b;
+      V vec1( P(0,0,0), random_sphere_point());
+      V vec2( P(0,0,0), random_sphere_point());
+      vec1*=0.9;
+      vec2*=0.9;
+      for(int i=0; i<N; ++i){
+        sp_a.push_back(random_sphere_point()+vec1);
+        sp_b.push_back(random_sphere_point()+vec2);
+      }
+      test(sp_a, sp_b, true);
+    }
+  }
+
+  void test_implicit_function(CGAL::Random &r, int N){
+    std::vector< Sphere<K> > spheres;
+    for(int i=0; i<N; ++i)
+      spheres.emplace_back(P(r.get_double(-2, 2),0,0),1);
+
+    for(int i=0; i<N; ++i)
+      for(int j=0; j<N; ++j)
+        assert(CGAL::Convex_hull_3::do_intersect(spheres[i], spheres[j])==(CGAL::abs(spheres[i].c.x()-spheres[j].c.x())<=2));
+  }
+
+  void full_test(CGAL::Random &r){
+    // test_degenerate();
+    // test_cube();
+    test_half_sphere();
+    test_random_tetrahedra(100, r);
+    // test_random_sphere(5000,20, r);
+  }
+
+};
+
+int main(int argc, char** argv)
+{
+  CGAL::Random r(argc==1?CGAL::get_default_random():std::stoi(argv[1]));
+  std::cout << "random seed = " << r.get_seed() << std::endl;
+
+  std::cout << std::setprecision(17);
+
+  // Test<Cartesian, CGAL::Surface_mesh<Cartesian::Point_3>>().full_test(r); //Too hard
+  Test<Epick, CGAL::Surface_mesh<Epick::Point_3> >().full_test(r);
+  Test<Epick, CGAL::Polyhedron_3<Epick> >().full_test(r);
+  Test<Epeck, CGAL::Surface_mesh<Epeck::Point_3> >().full_test(r);
+  return 0;
+}
diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib
index 7cfe999a4c6a..66206aeb58d5 100644
--- a/Documentation/doc/biblio/geom.bib
+++ b/Documentation/doc/biblio/geom.bib
@@ -152130,3 +152130,18 @@ @article{liu2025linequadrics
   url = {https://doi.org/10.1111/cgf.70184},
   year={2025}
 }
+
+@article{hornus2017convexintersection,
+  TITLE = {{Detecting the intersection of two convex shapes by searching on the 2-sphere}},
+  AUTHOR = {Hornus, Samuel},
+  URL = {https://inria.hal.science/hal-01522903},
+  NOTE = {Special issue: Proceedings of Solid and Physical Modeling (SPM 2017)},
+  JOURNAL = {{Computer-Aided Design}},
+  PUBLISHER = {{Elsevier}},
+  YEAR = {2017},
+  DOI = {10.1016/j.cad.2017.05.009},
+  KEYWORDS = {intersection detection ; GJK ; numerical robustness ; collision detection ; Minkowski sum ; Gauss map},
+  PDF = {https://inria.hal.science/hal-01522903v1/file/inter-spm.pdf},
+  HAL_ID = {hal-01522903},
+  HAL_VERSION = {v1},
+}
diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index 871b294b1120..87b9aa1712c1 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -38,6 +38,12 @@ Release date: July 2026
     parameter any longer. (This third optional parameter was introduced a few years ago, and now abandoned only for
     `do_intersect()`.)
 
+### [3D Convex hull](https://doc.cgal.org/6.2/Manual/packages.html#PkgConvexHull3)
+
+-   Added the functions `Convex_hull::do_intersect()` to the package `Convex_hull_3`, which enable testing the intersection of two convex hulls.
+-   Added the functions `extreme_point_3()` to the package `Convex_hull_3`, which return the farthest point of a convex hull in a given direction.
+-   Added the class `Convex_hull_hierarchy` to the package `Convex_hull_3`, which represents a convex hull and is optimized for use with the above functions.
+
 ## [Release 6.1](https://github.com/CGAL/cgal/releases/tag/v6.1)
 
 Release date: Sept 2025
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h b/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
index 42e84e436466..4230af91b5dc 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
@@ -862,6 +862,45 @@ const CGAL::Point_3<Kernel>& t);
 
 /// @}
 
+/// \defgroup compare_scalar_product_grp CGAL::compare_scalar_product()
+/// \ingroup kernel_global_function
+/// \sa `compare_distance_to_point_grp`
+/// \sa `compare_signed_distance_to_line_grp`
+/// \sa `compare_signed_distance_to_plane_grp`
+/// \sa `has_larger_distance_to_point_grp`
+/// \sa `has_larger_signed_distance_to_line_grp`
+/// \sa `has_larger_signed_distance_to_plane_grp`
+/// \sa `has_smaller_distance_to_point_grp`
+/// \sa `has_smaller_signed_distance_to_line_grp`
+/// \sa `has_smaller_signed_distance_to_plane_grp`
+/// @{
+
+/*!
+returns `CGAL::LARGER`
+iff the scalar product of `u` and
+`v` is larger than the scalar product `sp`, `CGAL::SMALLER`, iff it is smaller,
+and `CGAL::EQUAL` iff both are equal.
+*/
+template <typename Kernel>
+Comparison_result
+compare_scalar_product(const CGAL::Vector_3<Kernel> &u,
+                       const CGAL::Vector_3<Kernel> &v,
+                       const CGAL::FT<Kernel> &sp);
+
+/*!
+returns `CGAL::LARGER`
+iff the scalar product of `u` and
+`v` is larger than the scalar product of `w` and `x`, `CGAL::SMALLER`, iff it is smaller,
+and `CGAL::EQUAL` iff both are equal.
+*/
+template <typename Kernel>
+Comparison_result
+compare_scalar_product(const CGAL::Vector_3<Kernel> &u,
+                       const CGAL::Vector_3<Kernel> &v,
+                       const CGAL::Vector_3<Kernel> &w,
+                       const CGAL::Vector_3<Kernel> &x);
+
+/// @}
 
 
 /// \defgroup compare_slopes_grp CGAL::compare_slope()
diff --git a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
index c66bd0c922a3..c49114860e97 100644
--- a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
+++ b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
@@ -8988,6 +8988,45 @@ class LessSignedDistanceToPlane_3 {
 
 }; /* end Kernel::LessSignedDistanceToPlane_3 */
 
+/*!
+  \ingroup PkgKernel23ConceptsFunctionObjects
+  \cgalConcept
+
+  \cgalRefines{AdaptableTernaryFunction}
+
+  \sa `compare_scalar_product_grp`
+
+*/
+class CompareScalarProduct_3 {
+public:
+
+  /// \name Operations
+  /// A model of this concept must provide:
+  /// @{
+
+  /*!
+  returns `CGAL::LARGER` iff the scalar product of `u` and `v` is larger than the scalar product `sp`,
+  `CGAL::SMALLER`, iff it is smaller,
+  and `CGAL::EQUAL` iff both are equal.
+*/
+  Comparison_result operator()(const Kernel::Vector_3 &u,
+                               const Kernel::Vector_3 &v,
+                               const Kernel::FT &sp);
+
+  /*!
+    returns `CGAL::LARGER` iff the scalar product of `u` and `v` is larger than the scalar product of `w` and `x`,
+    `CGAL::SMALLER`, iff it is smaller,
+    and `CGAL::EQUAL` iff both are equal.
+  */
+  Comparison_result operator()(const Kernel::Vector_3 &u,
+                               const Kernel::Vector_3 &v,
+                               const Kernel::Vector_3 &w,
+                               const Kernel::Vector_3 &x);
+
+  /// @}
+
+}; /* end Kernel::CompareScalarProduct_3 */
+
 /*!
   \ingroup PkgKernel23ConceptsFunctionObjects
   \cgalConcept
diff --git a/Kernel_23/doc/Kernel_23/Concepts/Kernel.h b/Kernel_23/doc/Kernel_23/Concepts/Kernel.h
index e312763fd33f..d26da5321670 100644
--- a/Kernel_23/doc/Kernel_23/Concepts/Kernel.h
+++ b/Kernel_23/doc/Kernel_23/Concepts/Kernel.h
@@ -1465,6 +1465,11 @@ class Kernel {
   */
   typedef unspecified_type Less_signed_distance_to_plane_3;
 
+  /*!
+    a model of `Kernel::CompareScalarProduct_3`
+  */
+  typedef unspecified_type Compare_scalar_product_3;
+
   /*!
     a model of `Kernel::LessDistanceToPoint_3`
   */
diff --git a/Kernel_23/doc/Kernel_23/PackageDescription.txt b/Kernel_23/doc/Kernel_23/PackageDescription.txt
index f9b9d6991e7c..c216635cf08d 100644
--- a/Kernel_23/doc/Kernel_23/PackageDescription.txt
+++ b/Kernel_23/doc/Kernel_23/PackageDescription.txt
@@ -159,6 +159,7 @@
 - \link compare_lexicographically_linear_grp `CGAL::compare_lexicographically()` \endlink
 - \link compare_signed_distance_to_line_grp `CGAL::compare_signed_distance_to_line()` \endlink
 - \link compare_signed_distance_to_plane_grp `CGAL::compare_signed_distance_to_plane()` \endlink
+- \link compare_scalar_product_grp `CGAL::compare_scalar_product()` \endlink
 - \link compare_slopes_grp `CGAL::compare_slopes()` \endlink
 - `CGAL::compare_angle()`
 - \link compare_dihedral_angle_grp `CGAL::compare_dihedral_angle()` \endlink
diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h
index b71470a62625..72f8426f7b8b 100644
--- a/Kernel_23/include/CGAL/Kernel/function_objects.h
+++ b/Kernel_23/include/CGAL/Kernel/function_objects.h
@@ -187,6 +187,29 @@ namespace CommonKernelFunctors {
     { return assign(t, o); }
   };
 
+  template <typename K>
+  class Compare_scalar_product_3
+  {
+    typedef typename K::Comparison_result  Comparison_result;
+    typedef typename K::Vector_3           Vector_3;
+    typedef typename K::FT                 FT;
+  public:
+    Comparison_result
+    operator()(const Vector_3& u, const Vector_3& v, const FT& sp) const
+    {
+      typename K::Compute_scalar_product_3 scalar_product = K().compute_scalar_product_3_object();
+      return CGAL::compare(scalar_product(u,v), sp);
+    }
+
+    Comparison_result
+    operator()(const Vector_3& u, const Vector_3& v,
+               const Vector_3& w, const Vector_3& x) const
+    {
+      typename K::Compute_scalar_product_3 scalar_product = K().compute_scalar_product_3_object();
+      return CGAL::compare(scalar_product(u,v), scalar_product(w,x));
+    }
+  };
+
   template <typename K>
   class Compare_angle_3
   {
diff --git a/Kernel_23/include/CGAL/Kernel/global_functions_3.h b/Kernel_23/include/CGAL/Kernel/global_functions_3.h
index 4a2ef5d6b197..722b404857be 100644
--- a/Kernel_23/include/CGAL/Kernel/global_functions_3.h
+++ b/Kernel_23/include/CGAL/Kernel/global_functions_3.h
@@ -472,6 +472,27 @@ compare_lexicographically(const Point_3<K> &p,
     return internal::compare_lexicographically_xyz(p, q, K());
 }
 
+template < class K >
+inline
+typename K::Comparison_result
+compare_scalar_product(const Vector_3<K> &a,
+                       const Vector_3<K> &b,
+                       const typename K::FT &sp)
+{
+  return internal::compare_scalar_product(a, b, sp, K());
+}
+
+template < class K >
+inline
+typename K::Comparison_result
+compare_scalar_product(const Vector_3<K> &a,
+                       const Vector_3<K> &b,
+                       const Vector_3<K> &c,
+                       const Vector_3<K> &d)
+{
+  return internal::compare_scalar_product(a, b, c, d, K());
+}
+
 template < class K >
 inline
 typename K::Comparison_result
diff --git a/Kernel_23/include/CGAL/Kernel/global_functions_internal_3.h b/Kernel_23/include/CGAL/Kernel/global_functions_internal_3.h
index 9aa8e4a339ed..937e2d16b422 100644
--- a/Kernel_23/include/CGAL/Kernel/global_functions_internal_3.h
+++ b/Kernel_23/include/CGAL/Kernel/global_functions_internal_3.h
@@ -550,6 +550,30 @@ compare_lexicographically_xyz(const typename K::Point_3 &p,
   return k.compare_xyz_3_object()(p, q);
 }
 
+template < class K >
+inline
+typename K::Comparison_result
+compare_scalar_product(const typename K::Vector_3 &a,
+                       const typename K::Vector_3 &b,
+                       const typename K::FT &sp,
+                       const K &k)
+{
+  return k.compare_scalar_product_3_object()(a, b, sp);
+}
+
+template < class K >
+inline
+typename K::Comparison_result
+compare_scalar_product(const typename K::Vector_3 &a,
+                       const typename K::Vector_3 &b,
+                       const typename K::Vector_3 &c,
+                       const typename K::Vector_3 &d,
+                       const K &k)
+{
+  return k.compare_scalar_product_3_object()(a, b, c, d);
+}
+
+
 template < class K >
 inline
 typename K::Comparison_result
diff --git a/Kernel_23/include/CGAL/Kernel/interface_macros.h b/Kernel_23/include/CGAL/Kernel/interface_macros.h
index a2314aef1b79..5f1bb86aa35c 100644
--- a/Kernel_23/include/CGAL/Kernel/interface_macros.h
+++ b/Kernel_23/include/CGAL/Kernel/interface_macros.h
@@ -114,6 +114,8 @@ CGAL_Kernel_pred(Compare_power_distance_3,
                  compare_power_distance_3_object)
 CGAL_Kernel_pred(Compare_signed_distance_to_line_2,
                  compare_signed_distance_to_line_2_object)
+CGAL_Kernel_pred(Compare_scalar_product_3,
+                 compare_scalar_product_3_object)
 CGAL_Kernel_pred(Compare_slope_2,
                  compare_slope_2_object)
 CGAL_Kernel_pred(Compare_slope_3,
diff --git a/Kernel_23/test/Kernel_23/include/CGAL/_test_new_3.h b/Kernel_23/test/Kernel_23/include/CGAL/_test_new_3.h
index 4ce2d1975e6a..2dc8a8766af6 100644
--- a/Kernel_23/test/Kernel_23/include/CGAL/_test_new_3.h
+++ b/Kernel_23/test/Kernel_23/include/CGAL/_test_new_3.h
@@ -659,6 +659,16 @@ test_new_3(const R& rep)
   tmp34ab = CGAL::compare_distance(p1, p2, p3, p4);
   tmp34ab = CGAL::compare_distance(p1, p2, p3);
 
+typename R::Compare_scalar_product_3 compare_sp
+        = rep.compare_scalar_product_3_object();
+  tmp34ab = compare_sp(v3,v4,FT(1));
+  tmp34ab = compare_sp(v3,v4,v3,v4);
+  tmp34ab = compare_sp(v3,v4,v5,v3+v4);
+
+  tmp34ab = CGAL::compare_scalar_product(v3,v4,v3,v4);
+  tmp34ab = CGAL::compare_scalar_product(v3,v4,FT(1));
+  tmp34ab = CGAL::compare_scalar_product(v3,v4,v5,v3+v4);
+
   typename R::Compare_power_distance_3 compare_power_dist
         = rep.compare_power_distance_3_object();
   tmp34ab = compare_power_dist(p1, wp2, wp3);
diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
index 73c54858bf28..a9cfd2768775 100644
--- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
+++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
@@ -71,6 +71,7 @@ CGAL_add_named_parameter(vertex_feature_degree_t, vertex_feature_degree, vertex_
 
 // List of named parameters used in the package 'Polygon Mesh Processing'
 CGAL_add_named_parameter(geom_traits_t, geom_traits, geom_traits)
+CGAL_add_named_parameter(geom_traits_converter_t, geom_traits_converter, geom_traits_converter)
 CGAL_add_named_parameter(vertex_incident_patches_t, vertex_incident_patches, vertex_incident_patches_map)
 CGAL_add_named_parameter(density_control_factor_t, density_control_factor, density_control_factor)
 CGAL_add_named_parameter(use_delaunay_triangulation_t, use_delaunay_triangulation, use_delaunay_triangulation)
@@ -249,6 +250,9 @@ CGAL_add_named_parameter(number_of_relaxations_t, number_of_relaxations, number_
 // List of named parameters used in Optimal_bounding_box package
 CGAL_add_named_parameter(use_convex_hull_t, use_convex_hull, use_convex_hull)
 
+// List of named parameters used in Convex_hull_3 package
+CGAL_add_named_parameter(compute_convex_hull_t, compute_convex_hull, compute_convex_hull)
+
 // meshing parameters
 CGAL_add_named_parameter(boundary_subdivision_ratio_t, boundary_subdivision_ratio, boundary_subdivision_ratio)
 CGAL_add_named_parameter(subdivision_ratio_t, subdivision_ratio, subdivision_ratio)