Chapter 63
Spatial Sorting

Christophe Delage

Table of Contents

63.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.

63.2   Spatial Sorting

Cgal provides a small set of sorting algorithms, currently implemented only for 2D and 3D points, although it is easy to extend them to other objects through a traits mechanism.

Given an iterator range of points, the function hilbert_sort sorts them along the space-filling Hilbert curve. Combined with a std::random_shuffle, the function spatial_sort will sort points in a way keeping enough randomness to retain theoretical optimality for some algorithms¹, and close enough to a space-filling curve to speed-up algorithms.

The 2D and 3D triangulation classes of Cgal internally sort points as described above, when the points are passed as an iterator range to a constructor or the insert method.

63.3   Examples

In the following three examples you will see how to apply spatial sort to Cgal points, or on the indices into a container with points, or on a user defined point class.

63.3.1   Sorting Points

When you construct a triangulation from an iterator range of points, the points get spatially sorted internally.

If you have a need for inserting the points one by one you might sort them yourself.

File: examples/Spatial_sorting/example_delaunay_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/spatial_sort.h>

#include <iostream>
#include <iterator>
#include <algorithm>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> DT;

int main ()
{
    std::vector<K::Point_2> v;

    std::cout << "Reading..." << std::endl;
    std::istream_iterator<K::Point_2> begin(std::cin), end;
    std::copy(begin, end, std::back_inserter(v));

    // Comment the following three lines to get a massive speed down.
    std::cout << "Sorting..." << std::endl;
    std::random_shuffle(v.begin(), v.end());
    CGAL::spatial_sort(v.begin(), v.end());

    std::cout << "Delaunay..." << std::endl;
    DT dt;

    DT::Face_handle f;
    for (std::vector<K::Point_2>::const_iterator p = v.begin(); p != v.end(); ++p)
        f = dt.insert(*p, f)->face();

    std::cout << "Delaunay is_valid()..." << std::endl;
    dt.is_valid();

    std::cout << "Ok." << std::endl;

    return 0;
}

63.3.2   Sorting Indices

If you do not want to reorder your input, you can sort indices in your input instead.

File: examples/Spatial_sorting/sort_indices.cpp

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


typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_2 Point_2;
typedef std::vector<Point_2>::iterator Point_iterator;



template <typename Kernel, typename Iterator>
struct Sort_traits_2 {

  Kernel k;

  Sort_traits_2 (const Kernel &kernel = Kernel())
      : k (kernel)
  {}

  typedef Iterator Point_2;

  struct Less_x_2 {
    Kernel k;
    Less_x_2 (const Kernel &kernel = Kernel())
        : k (kernel)
    {}
    bool operator() (const Point_2 &p, const Point_2 &q) const
    {
      return k.less_x_2_object() (*p, *q);
    }
  };

  Less_x_2
  less_x_2_object() const
  {
    return Less_x_2(k);
  }

  struct Less_y_2 {
    Kernel k;
    Less_y_2 (const Kernel &kernel = Kernel())
        : k (kernel)
    {}
    bool operator() (const Point_2 &p, const Point_2 &q) const
    {
      return k.less_y_2_object() (*p, *q);
    }
  };


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


typedef CGAL::Hilbert_sort_2<Sort_traits_2<K, Point_iterator> > Hilbert_sort_2;
typedef CGAL::Multiscale_sort<Hilbert_sort_2> Spatial_sort_2;




int main ()
{
  Spatial_sort_2 sort_2;

  std::vector<Point_2> points;
  std::vector<Point_iterator> iterators;

  Point_2 p;
  while(std::cin >> p){
    points.push_back(p);
  }

  iterators.reserve(points.size());
  for(Point_iterator it = points.begin(); it != points.end(); ++it){
    iterators.push_back(it);
  }

  sort_2(iterators.begin(), iterators.end());

  for(std::vector<Point_iterator>::iterator i = iterators.begin();
          i != iterators.end(); i++)
  {
    std::cout << **i << std::endl;
  }

  return 0;
}

63.3.3   Sorting User Defined Points

If you want to sort points of your own point type, you only have to provide functors that compare the x and y coordinates of your points.

File: examples/Spatial_sorting/myPoint.cpp

#include <CGAL/spatial_sort.h>

struct MyPoint {
  double x,y;

  MyPoint()
    : x(0), y(0)
  {}

  MyPoint(double x, double y)
    : x(x), y(y)
  {}
};


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()
{
  MyPoint points[2];

  points[0] = MyPoint(78,12);
  points[1] = MyPoint(3,121);
  MySpatialSortingTraits sst;
  CGAL::spatial_sort(points, points+2, sst);
  std::cerr << "done" << std::endl;
  return 0;
}

Footnotes