Spatial_search: Duplicate point handling

Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp

@@ -22,6 +22,7 @@ struct Custom_traits_Hausdorff
     FT operator+(FT)const{return FT();}
     FT operator*(FT)const{return FT();}
     bool operator<=(FT)const{return false;}
+    bool operator==(FT)const{return false;}
     bool operator>=(FT)const{return false;}
     bool operator!=(FT)const{return false;}
     bool operator<(FT)const{return false;}

Spatial_searching/include/CGAL/Kd_tree.h

@@ -196,15 +196,15 @@ class Kd_tree {
     if (try_parallel_internal_node_creation (nh, c, c_low, tag))
       return;
 
-    if (c_low.size() > split.bucket_size())
+    if (c_low.size() > split.bucket_size() && !CGAL::is_zero(c_low.max_tight_spread()))
     {
       nh->lower_ch = new_internal_node();
       create_internal_node (nh->lower_ch, c_low, tag);
     }
     else
       nh->lower_ch = create_leaf_node(c_low);
 
-    if (c.size() > split.bucket_size())
+    if (c.size() > split.bucket_size() && !CGAL::is_zero(c.max_tight_spread()))
     {
       nh->upper_ch = new_internal_node();
       create_internal_node (nh->upper_ch, c, tag);
@@ -341,7 +341,7 @@ class Kd_tree {
 
     Point_container c(dim_, data.begin(), data.end(),traits_);
     bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box());
-    if (c.size() <= split.bucket_size()){
+    if (c.size() <= split.bucket_size() || CGAL::is_zero(c.max_tight_spread())){
       tree_root = create_leaf_node(c);
     }else {
        tree_root = new_internal_node();

Spatial_searching/include/CGAL/Point_container.h

@@ -419,7 +419,28 @@ class Point_container {
 
     typename Traits::Cartesian_const_iterator_d mpit = construct_it((*(*mid)));
     FT val1 = *(mpit+split_coord);
+
+    // Avoid using the low coord value as it results in an empty split
+    if (val1 == tbox.min_coord(split_coord)) {
+      iterator it = std::min_element(mid, end(), [=](const Point_d* a, const Point_d* b) -> bool {
+        FT a_c = *(construct_it(*a) + split_coord);
+        FT b_c = *(construct_it(*b) + split_coord);
+
+        if (a_c == val1)
+          return false;
+
+        if (b_c == val1)
+          return true;
+
+        return a_c < b_c;
+        });
+      return *(construct_it(**it) + split_coord);
+    }
+
     mid++;
+    if (mid == end())
+      return val1;
+
     mpit = construct_it((*(*mid)));
     FT val2 = *(mpit+split_coord);
     return (val1+val2)/FT(2);

Spatial_searching/include/CGAL/Splitters.h

@@ -19,6 +19,7 @@
 
 #include <CGAL/license/Spatial_searching.h>
 
+#include <CGAL/basic.h>
 #include <CGAL/Point_container.h>
 #include <CGAL/Plane_separator.h>
 
@@ -236,7 +237,9 @@ namespace CGAL {
     void
     operator() (Separator& sep, Container& c0, Container& c1) const
     {
-      sep = Separator(c0.max_span_coord(),FT(0));
+      if (!CGAL::is_zero(c0.max_spread()))
+        sep = Separator(c0.max_span_coord(),FT(0));
+      else sep = Separator(c0.max_tight_span_coord(), FT(0));
       sep.set_cutting_value(c0.median(sep.cutting_dimension()));
       c0.split(c1,sep,true);
     }