Alpha_kernel_d, its tests and Alpha complex to use it

6 files changed, 387 insertions, 49 deletions

diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index 2eecd50c..c9b937ed 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -7,10 +7,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
   target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY})
   add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp )
   target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY})
+  add_executable ( Weighted_alpha_complex_example_from_points Weighted_alpha_complex_from_points.cpp )
+  target_link_libraries(Weighted_alpha_complex_example_from_points  ${CGAL_LIBRARY})
   if (TBB_FOUND)
     target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES})
     target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES})
     target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES})
+    target_link_libraries(Weighted_alpha_complex_example_from_points ${TBB_LIBRARIES})
   endif()
 
   add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>)
diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp
new file mode 100644
index 00000000..19a04282
--- /dev/null
+++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp
@@ -0,0 +1,54 @@
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <CGAL/Epeck_d.h>
+
+#include <iostream>
+#include <vector>
+
+// Explicit dimension 2 Epeck_d kernel
+using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<2> >;
+using Bare_point = Kernel::Point_d;
+using Weighted_point = Kernel::Weighted_point_d;
+using Vector_of_points = std::vector<Weighted_point>;
+
+int main() {
+  // ----------------------------------------------------------------------------
+  // Init of a list of points
+  // ----------------------------------------------------------------------------
+  Vector_of_points points;
+  points.push_back(Weighted_point(Bare_point(1.0, 1.0) , 1.));
+  points.push_back(Weighted_point(Bare_point(7.0, 0.0) , 1.));
+  points.push_back(Weighted_point(Bare_point(4.0, 6.0) , 1.));
+  points.push_back(Weighted_point(Bare_point(9.0, 6.0) , 1.));
+  points.push_back(Weighted_point(Bare_point(0.0, 14.0), 1.));
+  points.push_back(Weighted_point(Bare_point(2.0, 19.0), 1.));
+  points.push_back(Weighted_point(Bare_point(9.0, 17.0), 1.));
+
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from the list of points
+  // ----------------------------------------------------------------------------
+  Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_weighted_points(points);
+
+  Gudhi::Simplex_tree<> simplex;
+  if (alpha_complex_from_weighted_points.create_complex(simplex)) {
+    // ----------------------------------------------------------------------------
+    // Display information about the alpha complex
+    // ----------------------------------------------------------------------------
+    std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() <<
+        " - " << simplex.num_simplices() << " simplices - " <<
+        simplex.num_vertices() << " vertices." << std::endl;
+
+    std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    for (auto f_simplex : simplex.filtration_simplex_range()) {
+      std::clog << "   ( ";
+      for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
+        std::clog << vertex << " ";
+      }
+      std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
+      std::clog << std::endl;
+    }
+  }
+  return 0;
+}
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index ba91998d..8e9fe773 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -12,6 +12,7 @@
 #ifndef ALPHA_COMPLEX_H_
 #define ALPHA_COMPLEX_H_
 
+#include <gudhi/Alpha_complex/Alpha_kernel_d.h>
 #include <gudhi/Debug_utils.h>
 // to construct Alpha_complex from a OFF file of points
 #include <gudhi/Points_off_io.h>
@@ -20,6 +21,7 @@
 #include <math.h>  // isnan, fmax
 
 #include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Regular_triangulation.h>  // aka. Weighted Delaunay triangulation
 #include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
 #include <CGAL/Epick_d.h>  // For FAST version
 #include <CGAL/Spatial_sort_traits_adapter_d.h>
@@ -91,49 +93,55 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val
  * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces).
  * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal &ge; 5.0.0.
  */
-template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>
+template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false>
 class Alpha_complex {
  public:
+  /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/
+  using Geom_traits = typename std::conditional<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>,
+                                                          Kernel>::type;
   // Add an int in TDS to save point index in the structure
-  typedef CGAL::Triangulation_data_structure<typename Kernel::Dimension,
-                              CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
-                              CGAL::Triangulation_full_cell<Kernel> > TDS;
-  /** \brief A Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
-  typedef CGAL::Delaunay_triangulation<Kernel, TDS> Delaunay_triangulation;
+  using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension,
+                                                 CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>,
+                                                 CGAL::Triangulation_full_cell<Geom_traits> >;
+  /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
+  using Triangulation = typename std::conditional<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
+                                                            CGAL::Delaunay_triangulation<Kernel, TDS>>::type;
+
+  using A_kernel_d = Alpha_kernel_d<Kernel, Weighted>;
+
+  // Numeric type of coordinates in the kernel
+  using FT = typename A_kernel_d::FT;
+
+  /** \brief If Weighted, the weighted point is cached (point + weight [= squared radius]),
+   * else a pair of point and squared radius is cached.
+  */
+  using Sphere = typename A_kernel_d::Sphere;
 
   /** \brief A point in Euclidean space.*/
-  typedef typename Kernel::Point_d Point_d;
-  /** \brief Geometric traits class that provides the geometric types and predicates needed by Delaunay
-   * triangulations.*/
-  typedef Kernel Geom_traits;
+  using Point_d = typename std::conditional<Weighted, typename A_kernel_d::Weighted_point_d,
+                                                      typename A_kernel_d::Point_d>::type;
 
  private:
-  typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
-  typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
-  typedef typename Kernel::Point_dimension_d        Point_Dimension;
-
   // Vertex_iterator type from CGAL.
-  typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
+  using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator;
 
   // size_type type from CGAL.
-  typedef typename Delaunay_triangulation::size_type size_type;
+  using size_type = typename Triangulation::size_type;
 
   // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
-  typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator;
-
-  // Numeric type of coordinates in the kernel
-  typedef typename Kernel::FT FT;
+  using Vector_vertex_iterator = typename std::vector< CGAL_vertex_iterator >;
 
  private:
   /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator.
    * Vertex handles are inserted sequentially, starting at 0.*/
   Vector_vertex_iterator vertex_handle_to_iterator_;
   /** \brief Pointer on the CGAL Delaunay triangulation.*/
-  Delaunay_triangulation* triangulation_;
+  Triangulation* triangulation_;
   /** \brief Kernel for triangulation_ functions access.*/
-  Kernel kernel_;
+  A_kernel_d kernel_;
+
   /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/
-  std::vector<std::pair<Point_d, FT>> cache_, old_cache_;
+  std::vector<Sphere> cache_, old_cache_;
 
  public:
   /** \brief Alpha_complex constructor from an OFF file name.
@@ -160,10 +168,10 @@ class Alpha_complex {
    *
    * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
    * 
-   * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d
+   * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d.
    * 
    * The type InputPointRange must be a range for which std::begin and
-   * std::end return input iterators on a Kernel::Point_d.
+   * std::end return input iterators on a Kernel::Point_d or Kernel::Weighted_point_d.
    */
   template<typename InputPointRange >
   Alpha_complex(const InputPointRange& points)
@@ -206,11 +214,8 @@ class Alpha_complex {
     auto last = std::end(points);
 
     if (first != last) {
-      // point_dimension function initialization
-      Point_Dimension point_dimension = kernel_.point_dimension_d_object();
-
-      // Delaunay triangulation is point dimension.
-      triangulation_ = new Delaunay_triangulation(point_dimension(*first));
+      // Delaunay triangulation init with point dimension.
+      triangulation_ = new Triangulation(kernel_.get_dimension(*first));
 
       std::vector<Point_d> point_cloud(first, last);
 
@@ -218,18 +223,22 @@ class Alpha_complex {
       std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
                                           boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
 
-      typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator,
-                                           CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map;
-      typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d;
+      using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator,
+                                                              CGAL::Identity_property_map<std::ptrdiff_t>>;
+      using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>;
 
       CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
 
-      typename Delaunay_triangulation::Full_cell_handle hint;
+      typename Triangulation::Full_cell_handle hint;
       for (auto index : indices) {
-        typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
-        // Save index value as data to retrieve it after insertion
-        pos->data() = index;
-        hint = pos->full_cell();
+        typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
+        if (pos != nullptr) {
+          // Save index value as data to retrieve it after insertion
+          pos->data() = index;
+          hint = pos->full_cell();
+        } else {
+          std::cout << "NULLPTR" << std::endl;
+        }
       }
       // --------------------------------------------------------------------------------------------
       // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
@@ -270,9 +279,7 @@ class Alpha_complex {
       v.clear();
       for (auto vertex : cplx.simplex_vertex_range(s))
         v.push_back(get_point_(vertex));
-      Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend());
-      FT r = kernel_.squared_distance_d_object()(c, v[0]);
-      cache_.emplace_back(std::move(c), std::move(r));
+      cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend()));
     }
     return cache_[k];
   }
@@ -282,13 +289,13 @@ class Alpha_complex {
   auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
     auto k = cplx.key(s);
     if(k!=cplx.null_key())
-      return old_cache_[k].second;
+      return kernel_.get_squared_radius(old_cache_[k]);
     // Using a transform_range is slower, currently.
     thread_local std::vector<Point_d> v;
     v.clear();
     for (auto vertex : cplx.simplex_vertex_range(s))
       v.push_back(get_point_(vertex));
-    return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend());
+    return kernel_.get_squared_radius(v.cbegin(), v.cend());
   }
 
  public:
@@ -322,9 +329,9 @@ class Alpha_complex {
                       bool exact = false,
                       bool default_filtration_value = false) {
     // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
-    typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
-    typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle;
-    typedef std::vector<Vertex_handle> Vector_vertex;
+    using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle;
+    using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle;
+    using Vector_vertex = std::vector<Vertex_handle>;
 
     if (triangulation_ == nullptr) {
       std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
@@ -416,8 +423,8 @@ class Alpha_complex {
   template <typename SimplicialComplexForAlpha, typename Simplex_handle>
   void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) {
     // From SimplicialComplexForAlpha type required to assign filtration values.
-    typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value;
-    typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
+    using Filtration_value = typename SimplicialComplexForAlpha::Filtration_value;
+    using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle;
 
     // ### Foreach Tau face of Sigma
     for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) {
@@ -450,7 +457,8 @@ class Alpha_complex {
         while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
         Vertex_handle extra = *longiter;
         auto const& cache=get_cache(complex, f_boundary);
-        bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second;
+        bool is_gab = kernel_.get_squared_distance(kernel_.get_circumcenter(cache), get_point_(extra)) >=
+                      kernel_.get_squared_radius(cache);
 #ifdef DEBUG_TRACES
         std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl;
 #endif  // DEBUG_TRACES
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h
new file mode 100644
index 00000000..87604ec4
--- /dev/null
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h
@@ -0,0 +1,137 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
+#define ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
+
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+#include <CGAL/Epick_d.h>  // For FAST version
+#include <CGAL/version.h>  // for CGAL_VERSION_NR
+
+#include <Eigen/src/Core/util/Macros.h>  // for EIGEN_VERSION_AT_LEAST
+
+#include <utility>  // for std::make_pair
+
+// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
+#if CGAL_VERSION_NR < 1041101000
+# error Alpha_complex is only available for CGAL >= 4.11
+#endif
+
+#if !EIGEN_VERSION_AT_LEAST(3,1,0)
+# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
+#endif
+
+namespace Gudhi {
+
+namespace alpha_complex {
+
+template < typename Kernel, bool Weighted = false >
+class Alpha_kernel_d {
+};
+
+// Unweighted Kernel_d version
+template < typename Kernel >
+class Alpha_kernel_d<Kernel, false> {
+ private:
+  // Kernel for functions access.
+  Kernel kernel_;
+ public:
+  // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h)
+  using Weighted_point_d = void;
+  using Point_d = typename Kernel::Point_d;
+  // Numeric type of coordinates in the kernel
+  using FT = typename Kernel::FT;
+  // Sphere is a pair of point and squared radius.
+  using Sphere = typename std::pair<Point_d, FT>;
+
+  int get_dimension(const Point_d& p0) const {
+    return kernel_.point_dimension_d_object()(p0);
+  }
+
+  template<class PointIterator>
+  Sphere get_sphere(PointIterator begin, PointIterator end) const {
+    Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
+    FT r = kernel_.squared_distance_d_object()(c, *begin);
+    return std::make_pair(std::move(c), std::move(r));
+  }
+
+  template<class PointIterator>
+  FT get_squared_radius(PointIterator begin, PointIterator end) const {
+    return kernel_.compute_squared_radius_d_object()(begin, end);
+  }
+
+  FT get_squared_radius(const Sphere& sph) const {
+    return sph.second;
+  }
+
+  template<class Point>
+  Point get_circumcenter(const Sphere& sph) const {
+    return sph.first;
+  }
+
+  template<class Point>
+  FT get_squared_distance(const Point& first, const Point& second) const {
+    return kernel_.squared_distance_d_object()(first, second);
+  }
+};
+
+// Weighted Kernel_d version
+template < typename Kernel >
+class Alpha_kernel_d<Kernel, true> {
+ private:
+  // Kernel for functions access.
+  Kernel kernel_;
+ public:
+  // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h)
+  using Point_d = void;
+  using Weighted_point_d = typename Kernel::Weighted_point_d;
+  using Bare_point_d = typename Kernel::Point_d;
+  // Numeric type of coordinates in the kernel
+  using FT = typename Kernel::FT;
+  // Sphere is a weighted point (point + weight [= squared radius]).
+  using Sphere = Weighted_point_d;
+
+  int get_dimension(const Weighted_point_d& p0) const {
+    return kernel_.point_dimension_d_object()(p0.point());
+  }
+
+  template<class PointIterator>
+  Sphere get_sphere(PointIterator begin, PointIterator end) const {
+    return kernel_.power_center_d_object()(begin, end);
+  }
+
+  template<class PointIterator>
+  FT get_squared_radius(PointIterator begin, PointIterator end) const {
+    Sphere sph = get_sphere(begin, end);
+    return sph.weight();
+  }
+
+  FT get_squared_radius(const Sphere& sph) const {
+    return sph.weight();
+  }
+
+  template<class Point>
+  Point get_circumcenter(const Sphere& sph) const {
+    return sph.point();
+  }
+
+  template<class Point>
+  FT get_squared_distance(const Point& first, const Point& second) const {
+    return kernel_.power_distance_d_object()(first, second);
+  }
+};
+
+}  // namespace alpha_complex
+
+namespace alphacomplex = alpha_complex;
+
+}  // namespace Gudhi
+
+#endif  // ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
\ No newline at end of file
diff --git a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
new file mode 100644
index 00000000..70b65aad
--- /dev/null
+++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
@@ -0,0 +1,132 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "alpha_kernel_d"
+#include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
+#include <CGAL/NT_converter.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>  // for std::pair
+
+#include <gudhi/Alpha_complex/Alpha_kernel_d.h>
+#include <gudhi/Unitary_tests_utils.h>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Exact_kernel_s;
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Inexact_kernel_s;
+// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
+
+typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_dimension, TestedKernel, list_of_kernel_variants) {
+  Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false> kernel;
+  std::vector<double> p0 {0., 1., 2., 3.};
+  typename TestedKernel::Point_d p0_d(p0.begin(), p0.end());
+
+  std::clog << "Dimension is " << kernel.get_dimension(p0_d) << std::endl;
+  BOOST_CHECK(kernel.get_dimension(p0_d) == 4);
+
+  Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true> w_kernel;
+  typename TestedKernel::Weighted_point_d w_p0_d(p0_d, 10.);
+
+  std::clog << "Dimension is " << w_kernel.get_dimension(w_p0_d) << std::endl;
+  BOOST_CHECK(w_kernel.get_dimension(w_p0_d) == 4);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kernel_variants) {
+  using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false>;
+  // Sphere: (x-1)² + (y-1)² + z² + t² = 1
+  // At least 5 points for a 3-sphere
+  std::vector<double> p0 {1., 0., 0., 0.};
+  std::vector<double> p1 {0., 1., 0., 0.};
+  std::vector<double> p2 {1., 1., 1., 0.};
+  std::vector<double> p3 {1., 1., 0., 1.};
+  std::vector<double> p4 {1., 1., -1., 0.};
+
+  using Point_d = typename Unweighted_kernel::Point_d;
+  std::vector<Point_d> unw_pts;
+  unw_pts.push_back(Point_d(p0.begin(), p0.end()));
+  unw_pts.push_back(Point_d(p1.begin(), p1.end()));
+  unw_pts.push_back(Point_d(p2.begin(), p2.end()));
+  unw_pts.push_back(Point_d(p3.begin(), p3.end()));
+  unw_pts.push_back(Point_d(p4.begin(), p4.end()));
+
+  Unweighted_kernel kernel;
+  auto unw_sphere = kernel.get_sphere(unw_pts.cbegin(), unw_pts.cend());
+
+  std::clog << "Center is " << unw_sphere.first << " - squared radius is " << unw_sphere.second << std::endl;
+
+  using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true>;
+
+  using Weighted_point_d = typename Weighted_kernel::Weighted_point_d;
+  using Bare_point_d = typename Weighted_kernel::Bare_point_d;
+  std::vector<Weighted_point_d> w_pts;
+  w_pts.push_back(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.));
+  w_pts.push_back(Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.));
+  w_pts.push_back(Weighted_point_d(Bare_point_d(p2.begin(), p2.end()), 0.));
+  w_pts.push_back(Weighted_point_d(Bare_point_d(p3.begin(), p3.end()), 0.));
+  w_pts.push_back(Weighted_point_d(Bare_point_d(p4.begin(), p4.end()), 0.));
+
+  Weighted_kernel w_kernel;
+  auto w_sphere = w_kernel.get_sphere(w_pts.cbegin(), w_pts.cend());
+
+  std::clog << "Center is " << w_sphere.point() << " - squared radius is " << w_sphere.weight() << std::endl;
+
+  CGAL::NT_converter<typename Weighted_kernel::FT, double> cast_to_double;
+  // The results shall be the same with weights = 0.
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(w_sphere.weight()));
+  BOOST_CHECK(unw_sphere.first == w_sphere.point());
+
+  auto unw_sq_rd = kernel.get_squared_radius(unw_pts.cbegin(), unw_pts.cend());
+  std::clog << "Squared radius is " << unw_sq_rd << std::endl;
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(unw_sq_rd));
+  auto w_sq_rd = w_kernel.get_squared_radius(w_pts.cbegin(), w_pts.cend());
+  std::clog << "Squared radius is " << w_sq_rd << std::endl;
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(w_sphere.weight()), cast_to_double(w_sq_rd));
+}
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_distance, TestedKernel, list_of_kernel_variants) {
+  using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false>;
+
+  std::vector<double> p0 {1., 0., 0., 0.};
+  std::vector<double> p1 {0., 1., 0., 0.};
+
+  using Point_d = typename Unweighted_kernel::Point_d;
+  Unweighted_kernel kernel;
+  auto dist_01 = kernel.get_squared_distance(Point_d(p0.begin(), p0.end()), Point_d(p1.begin(), p1.end()));
+  std::clog << "Distance is " << dist_01 << std::endl;
+
+  using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true>;
+
+  using Weighted_point_d = typename Weighted_kernel::Weighted_point_d;
+  using Bare_point_d = typename Weighted_kernel::Bare_point_d;
+  std::vector<Weighted_point_d> w_pts;
+
+  Weighted_kernel w_kernel;
+  auto w_dist_01 = w_kernel.get_squared_distance(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.),
+                                                 Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.));
+  std::clog << "Distance is " << w_dist_01 << std::endl;
+
+  CGAL::NT_converter<typename Weighted_kernel::FT, double> cast_to_double;
+  // The results shall be the same with weights = 0.
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(dist_01), cast_to_double(w_dist_01));
+}
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index fe4b23e4..bd61c0d8 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -10,13 +10,17 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
   target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY})
   add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp )
   target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY})
+  add_executable ( Alpha_kernel_d_test_unit Alpha_kernel_d_unit_test.cpp )
+  target_link_libraries(Alpha_kernel_d_test_unit ${CGAL_LIBRARY})
   if (TBB_FOUND)
     target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES})
     target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES})
+    target_link_libraries(Alpha_kernel_d_test_unit ${TBB_LIBRARIES})
   endif()
 
   gudhi_add_boost_test(Alpha_complex_test_unit)
   gudhi_add_boost_test(Delaunay_complex_test_unit)
+  gudhi_add_boost_test(Alpha_kernel_d_test_unit)
 
   add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp )
   target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY})