diff --git a/Frechet_distance/doc/Frechet_distance/PackageDescription.txt b/Frechet_distance/doc/Frechet_distance/PackageDescription.txt
index 5eff28044c5..3c3f829b1f8 100644
--- a/Frechet_distance/doc/Frechet_distance/PackageDescription.txt
+++ b/Frechet_distance/doc/Frechet_distance/PackageDescription.txt
@@ -5,8 +5,6 @@
 /// You can add some text here if you want a description at the group level
 /// \ingroup PkgFrechetDistanceRef
 
-/// \todo use Construct_point_d --> currently no concept in CGAL taking 2,3, or a range of coordinates
-
 /*!
 \addtogroup PkgFrechetDistanceRef
 
diff --git a/Frechet_distance/include/CGAL/Frechet_distance/internal/curve.h b/Frechet_distance/include/CGAL/Frechet_distance/internal/curve.h
index 50271283f0c..04d04c0cd1a 100644
--- a/Frechet_distance/include/CGAL/Frechet_distance/internal/curve.h
+++ b/Frechet_distance/include/CGAL/Frechet_distance/internal/curve.h
@@ -225,12 +225,27 @@ public:
 
     bool operator!=(Curve const& other) const { return !(*this == other); }
 
-    static IFT squared_distance(Point const& p, Point const& q, const Traits& traits)
+private:
+    // private as it will do nasty things if p is a temporary
+    static auto begin(const Point& p, const Traits& traits)
+    {
+      Construct_cartesian_const_iterator_d ccci = traits.construct_cartesian_const_iterator_d_object();
+      return ccci(p);
+    }
+
+    static auto begin(const std::array<FT,dimension>& p, const Traits&)
+    {
+      return p.begin();
+    }
+
+public:
+    template <class P1, class P2>
+    static IFT squared_distance(P1 const& p, P2 const& q, const Traits& traits)
     {
       IFT res(0);
-      Construct_cartesian_const_iterator_d ccci = traits.construct_cartesian_const_iterator_d_object();
 
-      auto itp = ccci(p), itq = ccci(q);
+      auto itp = begin(p, traits);
+      auto itq = begin(q, traits);
 
       for (int d=0;d<Traits::Dimension::value; ++d)
       {
@@ -239,7 +254,8 @@ public:
       }
       return res;
     }
-    static IFT distance(Point const& p, Point const& q, const Traits& traits)
+    template <class P1, class P2>
+    static IFT distance(P1 const& p, P2 const& q, const Traits& traits)
     {
       return sqrt(squared_distance(p,q, traits));
     }
@@ -268,8 +284,9 @@ public:
 
 
     template <class P>
-    Point interpolate_at(P const& pt) const
+    std::array<FT,dimension> interpolate_at(P const& pt) const
     {
+        std::array<distance_t, dimension> b_coords;
 //        assert(pt.getFraction() >= Lambda<Self>(0) &&
 //               pt.getFraction() <= Lambda<Self>(1));
         assert((
@@ -278,27 +295,25 @@ public:
 //            is_zero(pt.getFraction()))));
             pt.getFraction()==0)));
 
-//        if (is_zero(pt.getFraction())) {
-        if (pt.getFraction()==0) {
-            return point(pt.getPoint());
-        }
-
-        std::array<distance_t, dimension> b_coords;
-        auto fraction = pt.getFraction().approx;
-        distance_t one_m_f = distance_t(1) - fraction;
         Construct_cartesian_const_iterator_d ccci = traits_.construct_cartesian_const_iterator_d_object();
         auto itp = ccci(point(pt.getPoint()));
+
+//        if (is_zero(pt.getFraction())) {
+        if (pt.getFraction()==0) {
+          for (int d=0; d<dimension; ++d)
+            b_coords[d] = *itp++;
+          return b_coords;
+        }
+
+        auto fraction = pt.getFraction().approx;
+        distance_t one_m_f = distance_t(1) - fraction;
+
         auto itq = ccci(point(pt.getPoint()+1));
 
         for (int d=0; d<dimension; ++d)
           b_coords[d] = one_m_f * (*itp++) +  (*itq++) * fraction;
 
-        if constexpr (dimension==2)
-          return Point(b_coords[0], b_coords[1]);
-        else if constexpr (dimension==3)
-          return Point(b_coords[0], b_coords[1], b_coords[2]);
-        else
-          return Point(b_coords.begin(), b_coords.end());
+        return b_coords;
     }
 
     Point interpolate_at(PointID const& pt) const { return point(pt); }
diff --git a/Frechet_distance/test/Frechet_distance/Frechet_distance_minimal_traits_2.cpp b/Frechet_distance/test/Frechet_distance/Frechet_distance_minimal_traits_2.cpp
index a30a90b4b22..910e942bdf3 100644
--- a/Frechet_distance/test/Frechet_distance/Frechet_distance_minimal_traits_2.cpp
+++ b/Frechet_distance/test/Frechet_distance/Frechet_distance_minimal_traits_2.cpp
@@ -6,9 +6,7 @@ struct MinimalFrechetTraits {
   using FT = double;
 
   struct Point_d
-  {
-    Point_d(double, double) {}
-  };
+  {};
 
   struct Construct_bbox_d
   {