Chapter 65
Spatial Sorting

Christophe Delage and Olivier Devillers

Table of Contents

65.1   Introduction

Many geometric algorithms implemented in Cgal are incremental, and thus their speed is dependent on the order of insertion. This package provides sorting algorithms that may considerably improve running times of such algorithms.

The rationale is to sort objects along a space-filling curve so that two objects close geometrically will be close in the insertion order with high probability. That way, parts of a data structure that will be looked at during an insertion will probably have been looked at in a recent insertion, and thus probably will be in cache memory instead of main memory. As another side-effect, these sorting functions usually improve memory locality of the data structures produced by incremental algorithms, sometimes leading to speed ups in other algorithm using these data structures.

Some algorithms have a good complexity under randomized hypotheses which contradicts the idea of sorting the input using any sorting criterion. In such a case, it is possible to introduce just a bit of randomness to be able to combine the good randomized complexity and the good effects of locality [ACR03].

The predicates used by this package are comparisons between coordinates, thus there is no robustness issue involved here, for example to choose the arithmetic of the kernel.

65.2   Hilbert Sorting

In 2D, one can construct a space filling curve, that is a mapping f of [0,1] to the unit square [0,1]², such that f(0)=(0,0) and f(1)=(1,0) in the following way: the unit square is subdivided in four such that

f([0,(1)/(4)])=[0,(1)/(2)]², f([(1)/(4),(1)/(2)])=[0,(1)/(2)] × [(1)/(2),1], f([(1)/(2),(3)/(4)])=[(1)/(2),1]², and f([(3)/(4),1])=[(1)/(2),1] × [0,(1)/(2)].
f((1)/(4))=(0,(1)/(2))
f((1)/(2))=((1)/(2),(1)/(2)), and f((3)/(4))=(1,(1)/(2)).

Then each square is subdivided in the same way recursively. Figure 65.1 illustrates this process.

Now given a set of 2D points, they can be sorted in the order they have on such a space filling curve as illustrated in Figure 65.2 :

65.2.1   Examples

The code to use Hilbert sort is as simple as the following example:

File: examples/Spatial_sorting/hilbert.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/hilbert_sort.h>
#include <iostream>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2                                          Point;
typedef CGAL::Creator_uniform_2<double,Point>               Creator;

int main ()
{
  std::size_t size = 16;
  std::vector<Point> v; v.reserve(size);
  CGAL::points_on_square_grid_2(3.0, size,            // generate points
                                std::back_inserter(v), Creator());

  CGAL::hilbert_sort (v.begin(), v.end());            // sort

  for(std::size_t i=0; i<size; ++i)std::cout<<v[i]<<std::endl;//output
  return 0;
}

If instead of subdividing the square in a fixed way at its middle point, as above, we subdivide it by splitting at the median point (in x or y directions alternating), we construct a 2-d tree adapted to the point set. This tree can be visited in a similar manner and we get also a suitable ordering of the points (see Figure 65.3).

Cgal provides Hilbert sorting for points in 2D, 3D and higher dimensions, in the middle and the median policies.

The middle policy is easier to analyze, and is interesting in practice for well distributed set of points in small dimension (if the number of points is really smaller than 2^d). The median policy should be prefered for high dimension or if the point set distribution is not regular (or unknown). Since the median policy cannot be much worse than the middle policy, while the converse can happen, the median policy is the default behavior. Most theoretical results are using the middle policy [ACR03, But71, BG89, PB89].

This other example illustrates the use of the two different policies

File: examples/Spatial_sorting/hilbert_policies.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/hilbert_sort.h>
#include <iostream>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2                                          Point;

int main ()
{
  std::vector<Point> v; v.reserve(4);
  v.push_back( Point(0.0,0.0)) ;
  v.push_back( Point(1.0,1.0)) ;
  v.push_back( Point(0.1,0.1)) ;
  v.push_back( Point(0.2,0.8)) ;
  
  std::cout << "Hilbert sort (middle policy)." << std::endl;
  CGAL::hilbert_sort (v.begin(), v.end(), CGAL::Hilbert_sort_middle_policy());
  std::cout<<v[0]<<"; "<<v[1]<<"; "<<v[2]<<"; "<<v[3]<<"; "<<std::endl;
  std::cout << "Hilbert sort (median policy)." << std::endl;
  CGAL::hilbert_sort (v.begin(), v.end(), CGAL::Hilbert_sort_median_policy());
  std::cout<<v[0]<<"; "<<v[1]<<"; "<<v[2]<<"; "<<v[3]<<"; "<<std::endl;
  return 0;
}

65.3   Spatial Sorting

Hilbert sort cannot be used directly before feeding a randomized algorithm. Thus, the trick is to organize the point set in random buckets of increasing sizes, Hilbert sort being used only inside a bucket.

It has been proved, in the context of Delaunay triangulation, that such an order provides enough randomness to combine the advantages of a random order and a space filling curve order [ACR03].

Cgal provides spatial sorting for points in 2D, 3D and higher dimensions, with the middle and the median policies for Hilbert sort in the buckets.

65.3.1   Basic Example

The following example shows that, on particular input, spatial sort runs much faster than a bad order or than Hilbert sort (below results with release mode compilation on a 1.8GHz processor).

$ ./small_example_delaunay_2 
10000 points on a parabola
  Delaunay without spatial sort...         done in 6.33443 seconds.
  Delaunay with median hilbert sort...     done in 0.822975 seconds.
  Delaunay with median spatial sort...     done in 0.022415 seconds.

File: examples/Spatial_sorting/small_example_delaunay_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/spatial_sort.h>
#include <CGAL/Timer.h>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> DT;
void compute_delaunay(std::vector<K::Point_2>::iterator it, 
                        std::vector<K::Point_2>::iterator e){
    DT dt;
    DT::Face_handle hint;
    for( ;it!=e; ++it)  hint = dt.insert(*it, hint)->face();
}
int main ()
{   int size = 1000;
    std::vector<K::Point_2> v;
    v.reserve(size);
    CGAL::Timer cost;
    std::cout <<size<< " points on a parabola" << std::endl;
    for (int i=0; i< size; ++i) {
      double x= -size +i;
      v.push_back( K::Point_2( x, x*x ));
    }
    cost.reset();cost.start();
    std::cout << "  Delaunay without spatial sort... "<< std::flush;
    compute_delaunay(v.begin(),v.end());cost.stop();
    std::cout << "done in "<<cost.time()<<" seconds." << std::endl;
    cost.reset();cost.start();
    std::cout << "  Delaunay with Hilbert sort...    " << std::flush;
    CGAL::hilbert_sort(v.begin(),v.end());
    compute_delaunay(v.begin(),v.end());cost.stop();
    std::cout << "done in "<<cost.time()<<" seconds." << std::endl;
    cost.reset();cost.start();
    std::cout << "  Delaunay with spatial sort...    " << std::flush;
    CGAL::spatial_sort(v.begin(),v.end());
    compute_delaunay(v.begin(),v.end());cost.stop();
    std::cout << "done in "<<cost.time()<<" seconds." << std::endl;
    return 0;
}

65.3.2   Using Your Own Point Type

File: examples/Spatial_sorting/myPoint.cpp

#include <CGAL/spatial_sort.h>

struct MyPoint {
  double x,y;
  int color;
  MyPoint()
    : x(0), y(0),color(0)
  {}

  MyPoint(double x, double y, int color=0)
    : x(x), y(y), color(color)
  {}
};


struct MyLessX {

  bool operator()(const MyPoint& p, const MyPoint& q) const
  {
    return p.x < q.x;
  }

};

struct MyLessY {

  bool operator()(const MyPoint& p, const MyPoint& q) const
  {
    return p.y < q.y;
  }

};

struct MySpatialSortingTraits {

  typedef MyPoint Point_2;

  typedef MyLessX Less_x_2;
  typedef MyLessY Less_y_2;
  
  Less_x_2 less_x_2_object() const
  {
    return Less_x_2();
  }

  Less_y_2 less_y_2_object() const
  {
    return Less_y_2();
  }
};

int main()
{
  std::vector< MyPoint > points;

  points.push_back(MyPoint(14,12, 3));
  points.push_back(MyPoint(1,2  , 0));
  points.push_back(MyPoint(414,2, 5));
  points.push_back(MyPoint(4,21 , 1));
  points.push_back(MyPoint(7,74 , 2));
  points.push_back(MyPoint(74,4 , 4));  
  
  MySpatialSortingTraits sst;
  CGAL::spatial_sort(points.begin(), points.end(), sst);

  for (std::vector< MyPoint >::iterator it=points.begin();it!=points.end();++it)
    std::cout << it->color << " ";
  std::cout << "\n";  
  
  std::cerr << "done" << std::endl;
  return 0;
}

65.3.3   Sorting Arbitrary Types

The following examples illustrate the usage of these traits class adapters.

Sorting Using Pairs of Points and Integers

File: examples/Spatial_sorting/sp_sort_using_property_map_2.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/spatial_sort.h>
#include <CGAL/Spatial_sort_traits_adapter_2.h>
#include <vector>

typedef CGAL::Simple_cartesian<double>      Kernel;
typedef Kernel::Point_2                     Point_2;
typedef std::pair<Point_2,int>              Point_with_info;
typedef std::vector< Point_with_info >      Data_vector;

//property map
struct First_of_pair{
  //classical typedefs
  typedef Point_with_info key_type;
  typedef Point_2 value_type;
  typedef const Point_2& reference;
  typedef boost::readable_property_map_tag category;
};
//get function for property map
First_of_pair::reference
get(const First_of_pair&, const First_of_pair::key_type& k) {
  return k.first;
}

typedef CGAL::Spatial_sort_traits_adapter_2<Kernel,First_of_pair> Search_traits_2;

int main()
{
  Data_vector points;
  points.push_back(std::make_pair(Point_2(14,12) , 3));
  points.push_back(std::make_pair(Point_2(1,2)   , 0));
  points.push_back(std::make_pair(Point_2(414,2) , 5));
  points.push_back(std::make_pair(Point_2(4,21)  , 1));
  points.push_back(std::make_pair(Point_2(7,74)  , 2));
  points.push_back(std::make_pair(Point_2(74,4)  , 4));
  
  Search_traits_2 traits;
  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 << "done" << std::endl;
  
  return 0;
}

Sorting Using Indices of Points

File: examples/Spatial_sorting/sp_sort_using_property_map_3.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/spatial_sort.h>
#include <CGAL/Spatial_sort_traits_adapter_3.h>
#include <vector>
#include <boost/iterator/counting_iterator.hpp>

typedef CGAL::Simple_cartesian<double>                  Kernel;
typedef Kernel::Point_3                                 Point_3;
//using a pointer as a special property map type
typedef 
  CGAL::Spatial_sort_traits_adapter_3<Kernel,Point_3*>  Search_traits_3;

int main()
{
  std::vector<Point_3> points;
  points.push_back(Point_3(1,3,11));
  points.push_back(Point_3(14,34,46));
  points.push_back(Point_3(414,34,4));
  points.push_back(Point_3(4,2,56));
  points.push_back(Point_3(744,4154,43));
  points.push_back(Point_3(74,44,1));
  
  std::vector<std::ptrdiff_t> indices;
  indices.reserve(points.size());
  
  std::copy(boost::counting_iterator<std::ptrdiff_t>(0),
            boost::counting_iterator<std::ptrdiff_t>(points.size()),
            std::back_inserter(indices));
  
  CGAL::spatial_sort( indices.begin(),indices.end(),Search_traits_3(&(points[0])) );

  for (std::vector<std::ptrdiff_t>::iterator it=indices.begin();it!=indices.end();++it)
    std::cout << points[*it] << "\n";

  std::cout << "done" << std::endl;
  
  return 0;
}

Sorting Using Indices of Pairs of Points and Integers

File: examples/Spatial_sorting/sp_sort_using_property_map_d.cpp

#include <CGAL/Cartesian_d.h>
#include <CGAL/spatial_sort.h>
#include <CGAL/Spatial_sort_traits_adapter_d.h>
#include <boost/iterator/counting_iterator.hpp>
#include <vector>

typedef CGAL::Cartesian_d<double>           Kernel;
typedef Kernel::Point_d                     Point_d;
typedef std::pair<Point_d,int>              Point_with_info;
typedef std::vector< Point_with_info >      Data_vector;

//property map and get as friend
// to be allowed to use private member
class Vect_ppmap{
  const Data_vector& points;
public:
  //classical typedefs
  typedef Data_vector::size_type key_type;
  typedef Point_d value_type;
  typedef const value_type& reference;
  typedef boost::readable_property_map_tag category;

  Vect_ppmap(const Data_vector& points_):points(points_){}

  friend reference get(const Vect_ppmap& vmap, key_type i) {
    return vmap.points[i].first;
  }
};

typedef CGAL::Spatial_sort_traits_adapter_d<Kernel,Vect_ppmap>   Search_traits_d;

int main()
{
  double coords[] ={ 1.0, 1.0, 1.0, 1.0,
                     2.0, 2.0, 2.0, 2.0 };
  
  Data_vector points;
  points.push_back(std::make_pair(Point_d(4,coords  ,coords+4) , 1));
  points.push_back(std::make_pair(Point_d(4,coords+4,coords+8) , 2));

  std::vector<Vect_ppmap::key_type> indices;
  indices.reserve(points.size());  
  
  std::copy(
    boost::counting_iterator<Vect_ppmap::key_type>(0),
    boost::counting_iterator<Vect_ppmap::key_type>(points.size()),
    std::back_inserter(indices) );
  
  CGAL::spatial_sort( 
    indices.begin(),
    indices.end(),
    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 << " ";
  std::cout << std::endl;

  std::cout << "done" << std::endl;
  
  return 0;
}