diff --git a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_2.cpp b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_2.cpp
index 9d9145c5e96..129defe6842 100644
--- a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_2.cpp
+++ b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_2.cpp
@@ -24,11 +24,19 @@ int main()
   points.push_back(std::make_pair(Point_2(74,4)  , 4));
 
   Search_traits_2 traits;
+
+  std::cout << "Order using default policy (median)\n";
   CGAL::spatial_sort(points.begin(), points.end(), traits);
   for (Data_vector::iterator it=points.begin();it!=points.end();++it)
     std::cout << it->second << " ";
   std::cout << "\n";
 
+  std::cout << "Order using middle policy\n";
+  CGAL::spatial_sort(points.begin(), points.end(), traits, CGAL::Hilbert_sort_middle_policy());
+  for (Data_vector::iterator it=points.begin();it!=points.end();++it)
+    std::cout << it->second << " ";
+  std::cout << "\n";
+
   std::cout << "done" << std::endl;
 
   return 0;
diff --git a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_3.cpp b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_3.cpp
index 2e6f83847d4..7886931a66c 100644
--- a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_3.cpp
+++ b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_3.cpp
@@ -26,12 +26,22 @@ int main()
             boost::counting_iterator<std::size_t>(points.size()),
             std::back_inserter(indices));
 
+  std::cout << "Order using default policy (median)\n";
   CGAL::spatial_sort( indices.begin(),
                       indices.end(),
                       Search_traits_3(CGAL::make_property_map(points)) );
 
-  for (std::vector<std::size_t>::iterator it=indices.begin();it!=indices.end();++it)
-    std::cout << points[*it] << "\n";
+  for (std::size_t i : indices)
+    std::cout << points[i] << "\n";
+
+  std::cout << "Order using middle policy\n";
+  CGAL::spatial_sort( indices.begin(),
+                      indices.end(),
+                      Search_traits_3(CGAL::make_property_map(points)),
+                      CGAL::Hilbert_sort_middle_policy());
+
+  for (std::size_t i : indices)
+    std::cout << points[i] << "\n";
 
   std::cout << "done" << std::endl;
 
diff --git a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_d.cpp b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_d.cpp
index 382a366147b..19e9ce91883 100644
--- a/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_d.cpp
+++ b/Spatial_sorting/examples/Spatial_sorting/sp_sort_using_property_map_d.cpp
@@ -46,14 +46,27 @@ int main()
     boost::counting_iterator<Vect_ppmap::key_type>(points.size()),
     std::back_inserter(indices) );
 
+  std::cout << "Order using default policy (median)\n";
   CGAL::spatial_sort(
     indices.begin(),
     indices.end(),
-    Search_traits_d(Vect_ppmap(points)) );
+    Search_traits_d(Vect_ppmap(points))
+   );
 
-  std::vector<Vect_ppmap::key_type>::iterator it=indices.begin();
-  for (;it!=indices.end();++it)
-    std::cout << points[*it].second << " ";
+  for (auto i : indices)
+    std::cout << points[i].second << " ";
+  std::cout << std::endl;
+
+  std::cout << "Order using middle policy\n";
+  CGAL::spatial_sort(
+    indices.begin(),
+    indices.end(),
+    Search_traits_d(Vect_ppmap(points)),
+    CGAL::Hilbert_sort_middle_policy()
+   );
+
+  for (auto i : indices)
+    std::cout << points[i].second << " ";
   std::cout << std::endl;
 
   std::cout << "done" << std::endl;
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_median_2.h b/Spatial_sorting/include/CGAL/Hilbert_sort_median_2.h
index 23e87f7ac7b..722ccc3acbc 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_median_2.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_median_2.h
@@ -38,7 +38,7 @@ struct Hilbert_cmp_2<K,x,true>
 {
   typedef typename K::Point_2 Point;
   K k;
-  Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_2 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return Hilbert_cmp_2<K,x,false> (k) (q, p);
@@ -51,7 +51,7 @@ struct Hilbert_cmp_2<K,0,false>
 {
   typedef typename K::Point_2 Point;
   K k;
-  Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_2 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return k.less_x_2_object() (p, q);
@@ -64,7 +64,7 @@ struct Hilbert_cmp_2<K,1,false>
 {
   typedef typename K::Point_2 Point;
   K k;
-  Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_2 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return k.less_y_2_object() (p, q);
@@ -93,7 +93,7 @@ private:
   };
 
 public:
-  Hilbert_sort_median_2 (const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+  Hilbert_sort_median_2 (const Kernel &k, std::ptrdiff_t limit = 1)
     : _k(k), _limit (limit)
   {}
 
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_median_3.h b/Spatial_sorting/include/CGAL/Hilbert_sort_median_3.h
index 324a8bdb060..ebf71afed5a 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_median_3.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_median_3.h
@@ -35,7 +35,7 @@ struct Hilbert_cmp_3<K,x,true>
 {
   typedef typename K::Point_3 Point;
   K k;
-  Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_3 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return Hilbert_cmp_3<K,x,false> (k) (q, p);
@@ -48,7 +48,7 @@ struct Hilbert_cmp_3<K,0,false>
 {
   typedef typename K::Point_3 Point;
   K k;
-  Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_3 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return k.less_x_3_object() (p, q);
@@ -61,7 +61,7 @@ struct Hilbert_cmp_3<K,1,false>
 {
   typedef typename K::Point_3 Point;
   K k;
-  Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_3 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return k.less_y_3_object() (p, q);
@@ -74,7 +74,7 @@ struct Hilbert_cmp_3<K,2,false>
 {
   typedef typename K::Point_3 Point;
   K k;
-  Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+  Hilbert_cmp_3 (const K &_k) : k(_k) {}
   bool operator() (const Point &p, const Point &q) const
   {
     return k.less_z_3_object() (p, q);
@@ -109,7 +109,7 @@ private:
   };
 
 public:
-  Hilbert_sort_median_3 (const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+  Hilbert_sort_median_3 (const Kernel &k, std::ptrdiff_t limit = 1)
     : _k(k), _limit (limit)
   {}
 
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_median_d.h b/Spatial_sorting/include/CGAL/Hilbert_sort_median_d.h
index a88b3926006..44ddf9a7fa3 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_median_d.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_median_d.h
@@ -32,7 +32,7 @@ struct Hilbert_cmp_d
   K k;
   int axe;
   bool orient;
-  Hilbert_cmp_d (int a, bool o, const K &_k = K()) : k(_k), axe(a),  orient(o) {}
+  Hilbert_cmp_d (int a, bool o, const K &_k) : k(_k), axe(a),  orient(o) {}
 
   bool operator() (const Point &p, const Point &q) const
   {
@@ -64,7 +64,7 @@ private:
   };
 
 public:
-  Hilbert_sort_median_d(const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+  Hilbert_sort_median_d(const Kernel &k, std::ptrdiff_t limit = 1)
     : _k(k), _limit (limit)
   {}
 
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_2.h b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_2.h
index 32cf03cc435..1e1203f7d56 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_2.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_2.h
@@ -31,7 +31,7 @@ struct Fixed_hilbert_cmp_2<K,x,true>
   typedef typename K::Point_2 Point;
   K k;
   double value;
-  Fixed_hilbert_cmp_2 (double v, const K &_k = K()) : k(_k),value(v) {}
+  Fixed_hilbert_cmp_2 (double v, const K &_k) : k(_k),value(v) {}
   bool operator() (const Point &p) const
   {
     return ! Fixed_hilbert_cmp_2<K,x,false> (value, k) (p);
@@ -46,7 +46,7 @@ struct Fixed_hilbert_cmp_2<K,0,false>
   typedef typename K::Point_2 Point;
   K k;
   double value;
-  Fixed_hilbert_cmp_2 (double v, const K &_k = K()) : k(_k),value(v) {}
+  Fixed_hilbert_cmp_2 (double v, const K &_k) : k(_k),value(v) {}
   bool operator() (const Point &p) const
   {
     return to_double(k.compute_x_2_object()(p)) < value;
@@ -61,7 +61,7 @@ struct Fixed_hilbert_cmp_2<K,1,false>
   typedef typename K::Point_2 Point;
   K k;
   double value;
-  Fixed_hilbert_cmp_2 (double v, const K &_k = K()) : k(_k),value(v) {}
+  Fixed_hilbert_cmp_2 (double v, const K &_k) : k(_k),value(v) {}
   bool operator() (const Point &p) const
   {
     return to_double(k.compute_y_2_object()(p)) < value;
@@ -89,7 +89,7 @@ private:
   };
 
 public:
-  Hilbert_sort_middle_2 (const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+  Hilbert_sort_middle_2 (const Kernel &k, std::ptrdiff_t limit = 1)
     : _k(k), _limit (limit)
   {}
 
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_3.h b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_3.h
index 479cde7812d..a74e911515b 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_3.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_3.h
@@ -30,7 +30,7 @@ namespace internal {
         typedef typename K::Point_3 Point;
         K k;
         double value;
-        Fixed_hilbert_cmp_3 (double v, const K &_k = K()) : k(_k),value(v) {}
+        Fixed_hilbert_cmp_3 (double v, const K &_k) : k(_k),value(v) {}
         bool operator() (const Point &p) const
         {
           return ! Fixed_hilbert_cmp_3<K,x,false> (value,k) (p);
@@ -45,7 +45,7 @@ namespace internal {
         typedef typename K::Point_3 Point;
         K k;
         double value;
-        Fixed_hilbert_cmp_3 (double v, const K &_k = K()) : k(_k),value(v) {}
+        Fixed_hilbert_cmp_3 (double v, const K &_k) : k(_k),value(v) {}
         bool operator() (const Point &p) const
         {
           return to_double(k.compute_x_3_object()(p)) < value;
@@ -60,7 +60,7 @@ namespace internal {
         typedef typename K::Point_3 Point;
         K k;
         double value;
-        Fixed_hilbert_cmp_3 (double v, const K &_k = K()) : k(_k),value(v) {}
+        Fixed_hilbert_cmp_3 (double v, const K &_k) : k(_k),value(v) {}
         bool operator() (const Point &p) const
         {
           return to_double(k.compute_y_3_object()(p)) < value;
@@ -75,7 +75,7 @@ namespace internal {
         typedef typename K::Point_3 Point;
         K k;
         double value;
-        Fixed_hilbert_cmp_3 (double v, const K &_k = K()) : k(_k),value(v) {}
+        Fixed_hilbert_cmp_3 (double v, const K &_k) : k(_k),value(v) {}
         bool operator() (const Point &p) const
         {
           return to_double(k.compute_z_3_object()(p)) < value ;
@@ -98,7 +98,7 @@ private:
       { Cmp (double v,const Kernel &k) : internal::Fixed_hilbert_cmp_3<Kernel,x,up> (v,k) {} };
 
 public:
-    Hilbert_sort_middle_3 (const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+    Hilbert_sort_middle_3 (const Kernel &k, std::ptrdiff_t limit = 1)
         : _k(k), _limit (limit)
     {}
 
@@ -154,26 +154,25 @@ public:
     template <class RandomAccessIterator>
     void operator() (RandomAccessIterator begin, RandomAccessIterator end) const
     {
-      K k;
-      double xmin=to_double(k.compute_x_3_object()(*begin)),
-             ymin=to_double(k.compute_y_3_object()(*begin)),
-             zmin=to_double(k.compute_z_3_object()(*begin)),
+      double xmin=to_double(_k.compute_x_3_object()(*begin)),
+             ymin=to_double(_k.compute_y_3_object()(*begin)),
+             zmin=to_double(_k.compute_z_3_object()(*begin)),
                xmax=xmin,
              ymax=ymin,
              zmax=zmin;
       for(RandomAccessIterator it=begin+1; it<end; ++it){
-        if ( to_double(k.compute_x_3_object()(*it)) < xmin)
-          xmin = to_double(k.compute_x_3_object()(*it));
-        if ( to_double(k.compute_y_3_object()(*it)) < ymin)
-          ymin = to_double(k.compute_y_3_object()(*it));
-        if ( to_double(k.compute_z_3_object()(*it)) < zmin)
-          zmin = to_double(k.compute_z_3_object()(*it));
-        if ( to_double(k.compute_x_3_object()(*it)) > xmax)
-          xmax = to_double(k.compute_x_3_object()(*it));
-        if ( to_double(k.compute_y_3_object()(*it)) > ymax)
-          ymax = to_double(k.compute_y_3_object()(*it));
-        if ( to_double(k.compute_z_3_object()(*it)) > zmax)
-          zmax = to_double(k.compute_z_3_object()(*it));
+        if ( to_double(_k.compute_x_3_object()(*it)) < xmin)
+          xmin = to_double(_k.compute_x_3_object()(*it));
+        if ( to_double(_k.compute_y_3_object()(*it)) < ymin)
+          ymin = to_double(_k.compute_y_3_object()(*it));
+        if ( to_double(_k.compute_z_3_object()(*it)) < zmin)
+          zmin = to_double(_k.compute_z_3_object()(*it));
+        if ( to_double(_k.compute_x_3_object()(*it)) > xmax)
+          xmax = to_double(_k.compute_x_3_object()(*it));
+        if ( to_double(_k.compute_y_3_object()(*it)) > ymax)
+          ymax = to_double(_k.compute_y_3_object()(*it));
+        if ( to_double(_k.compute_z_3_object()(*it)) > zmax)
+          zmax = to_double(_k.compute_z_3_object()(*it));
       }
 
       sort <0, false, false, false> (begin, end, xmin,ymin,zmin,xmax,ymax,zmax);
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_d.h b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_d.h
index f802f6b2a1b..6d45051100e 100644
--- a/Spatial_sorting/include/CGAL/Hilbert_sort_middle_d.h
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_middle_d.h
@@ -32,7 +32,7 @@ namespace internal {
         int axe;
         bool orient;
         double value;
-        Fixed_hilbert_cmp_d (int a, bool o, double v, const K &_k = K())
+        Fixed_hilbert_cmp_d (int a, bool o, double v, const K &_k)
           :  k(_k), axe(a), orient(o), value(v) {}
         bool operator() (const Point &p) const
         {
@@ -64,7 +64,7 @@ private:
         : internal::Fixed_hilbert_cmp_d<Kernel> (a,dir,v,k) {} };
 
 public:
-    Hilbert_sort_middle_d (const Kernel &k = Kernel(), std::ptrdiff_t limit = 1)
+    Hilbert_sort_middle_d (const Kernel &k, std::ptrdiff_t limit = 1)
         : _k(k), _limit (limit)
     {}
 
diff --git a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_2.h b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_2.h
index bd1e9069693..e83efb6ceaf 100644
--- a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_2.h
+++ b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_2.h
@@ -37,6 +37,28 @@ public:
   typedef typename boost::property_traits<PointPropertyMap>::key_type Point_2;
   typedef typename boost::call_traits<Point_2>::param_type Arg_type;
 
+  struct Compute_x_2
+    : public Base_traits::Compute_x_2
+  {
+    Compute_x_2(const PointPropertyMap& ppmap, const typename Base_traits::Compute_x_2& base):
+      Base_traits::Compute_x_2(base), ppmap_(ppmap){}
+    const PointPropertyMap& ppmap_;
+    typename Gt::FT operator()(Arg_type p) const {
+      return static_cast<const typename Base_traits::Compute_x_2*>(this)->operator()(get(ppmap_,p));
+    }
+  };
+
+  struct Compute_y_2
+    : public Base_traits::Compute_y_2
+  {
+    Compute_y_2(const PointPropertyMap& ppmap, const typename Base_traits::Compute_y_2& base):
+      Base_traits::Compute_y_2(base), ppmap_(ppmap){}
+    const PointPropertyMap& ppmap_;
+    typename Gt::FT operator()(Arg_type p) const {
+      return static_cast<const typename Base_traits::Compute_y_2*>(this)->operator()(get(ppmap_,p));
+    }
+  };
+
   struct Less_x_2
     : public Base_traits::Less_x_2
   {
@@ -61,6 +83,8 @@ public:
 
   Less_x_2 less_x_2_object () const {return Less_x_2(ppmap_, static_cast<const Gt*>(this)->less_x_2_object() );}
   Less_y_2 less_y_2_object () const {return Less_y_2(ppmap_, static_cast<const Gt*>(this)->less_y_2_object() );}
+  Compute_x_2 compute_x_2_object () const {return Compute_x_2(ppmap_, static_cast<const Gt*>(this)->compute_x_2_object() );}
+  Compute_y_2 compute_y_2_object () const {return Compute_y_2(ppmap_, static_cast<const Gt*>(this)->compute_y_2_object() );}
 
   const PointPropertyMap& point_property_map() const {return ppmap_;}
 };
diff --git a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_3.h b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_3.h
index 778c392b09e..cc304d1965b 100644
--- a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_3.h
+++ b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_3.h
@@ -42,7 +42,7 @@ public:
     Compute_x_3(const PointPropertyMap& ppmap, const typename Base_traits::Compute_x_3& base):
       Base_traits::Compute_x_3(base), ppmap_(ppmap){}
     const PointPropertyMap& ppmap_;
-    bool operator()(Arg_type p) const {
+    typename Gt::FT operator()(Arg_type p) const {
       return static_cast<const typename Base_traits::Compute_x_3*>(this)->operator()(get(ppmap_,p));
     }
   };
@@ -53,7 +53,7 @@ public:
     Compute_y_3(const PointPropertyMap& ppmap, const typename Base_traits::Compute_y_3& base):
       Base_traits::Compute_y_3(base), ppmap_(ppmap){}
     const PointPropertyMap& ppmap_;
-    bool operator()(Arg_type p) const {
+    typename Gt::FT operator()(Arg_type p) const {
       return static_cast<const typename Base_traits::Compute_y_3*>(this)->operator()(get(ppmap_,p));
     }
   };
@@ -64,7 +64,7 @@ public:
     Compute_z_3(const PointPropertyMap& ppmap, const typename Base_traits::Compute_z_3& base):
       Base_traits::Compute_z_3(base), ppmap_(ppmap){}
     const PointPropertyMap& ppmap_;
-    bool operator()(Arg_type p) const {
+    typename Gt::FT operator()(Arg_type p) const {
       return static_cast<const typename Base_traits::Compute_z_3*>(this)->operator()(get(ppmap_,p));
     }
   };
diff --git a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_d.h b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_d.h
index 0a6313c244f..42be31f946f 100644
--- a/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_d.h
+++ b/Spatial_sorting/include/CGAL/Spatial_sort_traits_adapter_d.h
@@ -64,7 +64,7 @@ public:
     Compute_coordinate_d(const PointPropertyMap& ppmap, const typename Base_traits::Compute_coordinate_d& base):
       Base_traits::Compute_coordinate_d(base), ppmap_(ppmap){}
     const PointPropertyMap& ppmap_;
-    bool operator()(Arg_type p, int i) const {
+    typename Gt::FT operator()(Arg_type p, int i) const {
       return static_cast<const typename Base_traits::Compute_coordinate_d*>(this)->operator()(get(ppmap_,p), i);
     }
   };