diff --git a/.gitattributes b/.gitattributes
index 75f7a9761ab..2534e8a2ba7 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1304,6 +1304,14 @@ Snap_rounding_2/doc_tex/Snap_rounding_2/sr1.ps -text
 Spatial_searching/demo/Spatial_searching/spatial_searching.vcproj -text
 Spatial_searching/doc_tex/Spatial_searching/Fig1.gif -text svneol=unset#unset
 Spatial_searching/doc_tex/Spatial_searching/Fig1.ps -text
+Spatial_sorting/description.txt -text
+Spatial_sorting/include/CGAL/BRIO_sort.h -text
+Spatial_sorting/include/CGAL/Hilbert_sort_2.h -text
+Spatial_sorting/include/CGAL/Hilbert_sort_3.h -text
+Spatial_sorting/include/CGAL/Hilbert_sort_base.h -text
+Spatial_sorting/maintainer -text
+Spatial_sorting/test/Spatial_sorting/test_brio.C -text
+Spatial_sorting/test/Spatial_sorting/test_hilbert.C -text
 Straight_skeleton_2/demo/Straight_skeleton_2/data/alley_0.poly -text
 Straight_skeleton_2/demo/Straight_skeleton_2/data/alley_1.poly -text
 Straight_skeleton_2/demo/Straight_skeleton_2/data/alley_2.poly -text
diff --git a/Spatial_sorting/description.txt b/Spatial_sorting/description.txt
new file mode 100644
index 00000000000..7d453557198
--- /dev/null
+++ b/Spatial_sorting/description.txt
@@ -0,0 +1,2 @@
+Package Spatial_sorting:
+Clever sorting of points to improve incremental algorithms (typically Delaunay).
diff --git a/Spatial_sorting/include/CGAL/BRIO_sort.h b/Spatial_sorting/include/CGAL/BRIO_sort.h
new file mode 100644
index 00000000000..1500ec85348
--- /dev/null
+++ b/Spatial_sorting/include/CGAL/BRIO_sort.h
@@ -0,0 +1,72 @@
+#ifndef CGAL_BRIO_SORT_H
+#define CGAL_BRIO_SORT_H
+
+#include <CGAL/basic.h>
+
+#include <CGAL/Hilbert_sort_2.h>
+#include <CGAL/Hilbert_sort_3.h>
+
+CGAL_BEGIN_NAMESPACE
+
+template <class Sort, class RandomAccessIterator>
+class BRIO_sort
+{
+    Sort sort;
+    int threshold;
+    int ratio;
+public:
+    BRIO_sort (int _ratio = 8, int _threshold = 64, const Sort &_sort = Sort())
+        : sort (_sort), threshold (_threshold), ratio (_ratio) 
+    {}
+    
+    void operator() (RandomAccessIterator begin, RandomAccessIterator end) const
+    {
+        RandomAccessIterator middle = begin;
+        if (end - begin >= threshold) {
+            middle = begin + (end - begin) / ratio;
+            this->operator() (begin, middle);
+        }
+        sort (middle, end);
+    }
+};
+
+template <class RandomAccessIterator, class Kernel>
+void brio_sort_2 (RandomAccessIterator begin, RandomAccessIterator end, Kernel k)
+{
+    typedef Hilbert_sort_2<Kernel, RandomAccessIterator> Sort;
+    (BRIO_sort<Sort,RandomAccessIterator> (4, 16, Sort (k))) (begin, end);
+}
+
+template <class RandomAccessIterator>
+void brio_sort_2 (RandomAccessIterator begin, RandomAccessIterator end)
+{
+    typedef std::iterator_traits<RandomAccessIterator> ITraits;
+    typedef typename ITraits::value_type               value_type;
+    typedef CGAL::Kernel_traits<value_type>            KTraits;
+    typedef typename KTraits::Kernel                   Kernel;
+
+    brio_sort_2 (begin, end, Kernel ());
+}
+
+template <class RandomAccessIterator, class Kernel>
+void brio_sort_3 (RandomAccessIterator begin, RandomAccessIterator end, Kernel k)
+{
+    typedef Hilbert_sort_3<Kernel, RandomAccessIterator> Sort;
+    (BRIO_sort<Sort,RandomAccessIterator> (8, 64, Sort (k))) (begin, end);
+}
+
+template <class RandomAccessIterator>
+void brio_sort_3 (RandomAccessIterator begin, RandomAccessIterator end)
+{
+    typedef std::iterator_traits<RandomAccessIterator> ITraits;
+    typedef typename ITraits::value_type               value_type;
+    typedef CGAL::Kernel_traits<value_type>            KTraits;
+    typedef typename KTraits::Kernel                   Kernel;
+
+    brio_sort_3 (begin, end, Kernel ());
+}
+
+CGAL_END_NAMESPACE
+
+
+#endif//CGAL_BRIO_SORT_H
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_2.h b/Spatial_sorting/include/CGAL/Hilbert_sort_2.h
new file mode 100644
index 00000000000..b6bc598e2a4
--- /dev/null
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_2.h
@@ -0,0 +1,120 @@
+#ifndef CGAL_HILBERT_SORT_2_H
+#define CGAL_HILBERT_SORT_2_H
+
+#include <CGAL/basic.h>
+
+#include <CGAL/Hilbert_sort_base.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+    template <class K, int x, bool up> struct Hilbert_cmp_2;
+
+
+    template <class K, int x>
+    struct Hilbert_cmp_2<K,x,true>
+        : public std::binary_function<typename K::Point_2,
+                                      typename K::Point_2, bool>
+    {
+        typedef typename K::Point_2 Point;
+        K k;
+        Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        { 
+            return Hilbert_cmp_2<K,x,false> (k) (q, p);
+        }
+    };
+    
+    template <class K>
+    struct Hilbert_cmp_2<K,0,false>
+        : public std::binary_function<typename K::Point_2,
+                                      typename K::Point_2, bool>
+    {
+        typedef typename K::Point_2 Point;
+        K k;
+        Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        { 
+            return k.less_x_2_object() (p, q);
+        }
+    };
+    
+    template <class K>
+    struct Hilbert_cmp_2<K,1,false>
+        : public std::binary_function<typename K::Point_2,
+                                      typename K::Point_2, bool>
+    {
+        typedef typename K::Point_2 Point;
+        K k;
+        Hilbert_cmp_2 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        { 
+            return k.less_y_2_object() (p, q);
+        }
+    };
+};
+
+template <class K, class RandomAccessIterator>
+class Hilbert_sort_2
+{
+public:
+    typedef RandomAccessIterator Iterator;
+    typedef K Kernel;
+    typedef typename Kernel::Point_2 Point;
+    
+private:
+    Kernel k;
+
+    template <int x, bool up> struct Cmp : public CGALi::Hilbert_cmp_2<Kernel,x,up>
+    { Cmp (const Kernel &k) : CGALi::Hilbert_cmp_2<Kernel,x,up> (k) {} };
+
+public:
+    Hilbert_sort_2 (const Kernel &_k = Kernel()) : k(_k) {}
+
+    template <int x, bool upx, bool upy>
+    void sort (Iterator begin, Iterator end) const
+    {
+        const int y = (x + 1) % 2;
+        if (end - begin <= 8) return;
+
+        Iterator m0 = begin, m4 = end;
+
+        Iterator m2 = CGALi::hilbert_split (m0, m4, Cmp< x,  upx> (k));
+        Iterator m1 = CGALi::hilbert_split (m0, m2, Cmp< y,  upy> (k));
+        Iterator m3 = CGALi::hilbert_split (m2, m4, Cmp< y, !upy> (k));
+
+        if (end - begin <= 8*4) return;
+
+        sort<y, upy, upx> (m0, m1);
+        sort<x, upx, upy> (m1, m2);
+        sort<x, upx, upy> (m2, m3);
+        sort<y,!upy,!upx> (m3, m4);
+    }
+
+    void operator() (Iterator begin, Iterator end) const
+    {
+        sort <0, false, false> (begin, end);
+    }
+};
+
+template <class RandomAccessIterator, class Kernel>
+void hilbert_sort_2 (RandomAccessIterator begin, RandomAccessIterator end,
+                     Kernel k)
+{
+    (Hilbert_sort_2<Kernel,RandomAccessIterator> (k)) (begin, end);
+}
+
+template <class RandomAccessIterator>
+void hilbert_sort_2 (RandomAccessIterator begin, RandomAccessIterator end)
+{
+    typedef std::iterator_traits<RandomAccessIterator> ITraits;
+    typedef typename ITraits::value_type               value_type;
+    typedef CGAL::Kernel_traits<value_type>            KTraits;
+    typedef typename KTraits::Kernel                   Kernel;
+
+    hilbert_sort_2 (begin, end, Kernel());
+}
+
+CGAL_END_NAMESPACE
+
+#endif//CGAL_HILBERT_SORT_2_H
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_3.h b/Spatial_sorting/include/CGAL/Hilbert_sort_3.h
new file mode 100644
index 00000000000..d673d494b2b
--- /dev/null
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_3.h
@@ -0,0 +1,143 @@
+#ifndef CGAL_HILBERT_SORT_3_H
+#define CGAL_HILBERT_SORT_3_H
+
+#include <CGAL/basic.h>
+
+#include <CGAL/Hilbert_sort_base.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+    template <class K, int x, bool up> struct Hilbert_cmp_3;
+
+
+    template <class K, int x>
+    struct Hilbert_cmp_3<K,x,true>
+        : public std::binary_function<typename K::Point_3,
+                                      typename K::Point_3, bool>
+    {
+        typedef typename K::Point_3 Point;
+        K k;
+        Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        {
+            return Hilbert_cmp_3<K,x,false> (k) (q, p);
+        }
+    };
+
+    template <class K>
+    struct Hilbert_cmp_3<K,0,false>
+        : public std::binary_function<typename K::Point_3,
+                                      typename K::Point_3, bool>
+    {
+        typedef typename K::Point_3 Point;
+        K k;
+        Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        {
+            return k.less_x_3_object() (p, q);
+        }
+    };
+
+    template <class K>
+    struct Hilbert_cmp_3<K,1,false>
+        : public std::binary_function<typename K::Point_3,
+                                      typename K::Point_3, bool>
+    {
+        typedef typename K::Point_3 Point;
+        K k;
+        Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        {
+            return k.less_y_3_object() (p, q);
+        }
+    };
+
+    template <class K>
+    struct Hilbert_cmp_3<K,2,false>
+        : public std::binary_function<typename K::Point_3,
+                                      typename K::Point_3, bool>
+    {
+        typedef typename K::Point_3 Point;
+        K k;
+        Hilbert_cmp_3 (const K &_k = K()) : k(_k) {}
+        bool operator() (const Point &p, const Point &q) const
+        {
+            return k.less_z_3_object() (p, q);
+        }
+    };
+};
+
+
+template <class K, class RandomAccessIterator>
+class Hilbert_sort_3
+{
+public:
+    typedef RandomAccessIterator Iterator;
+    typedef K Kernel;
+    typedef typename Kernel::Point_3 Point;
+
+private:
+    Kernel k;
+
+    template <int x, bool up> struct Cmp : public CGALi::Hilbert_cmp_3<Kernel,x,up>
+    { Cmp (const Kernel &k) : CGALi::Hilbert_cmp_3<Kernel,x,up> (k) {} };
+
+public:
+    Hilbert_sort_3 (const Kernel &_k = Kernel()) : k(_k) {}
+
+    template <int x, bool upx, bool upy, bool upz>
+    void sort (Iterator begin, Iterator end) const
+    {
+        const int y = (x + 1) % 3, z = (x + 2) % 3;
+        if (end - begin <= 8) return;
+
+        Iterator m0 = begin, m8 = end;
+
+        Iterator m4 = CGALi::hilbert_split (m0, m8, Cmp< x,  upx> (k));
+        Iterator m2 = CGALi::hilbert_split (m0, m4, Cmp< y,  upy> (k));
+        Iterator m1 = CGALi::hilbert_split (m0, m2, Cmp< z,  upz> (k));
+        Iterator m3 = CGALi::hilbert_split (m2, m4, Cmp< z, !upz> (k));
+        Iterator m6 = CGALi::hilbert_split (m4, m8, Cmp< y, !upy> (k));
+        Iterator m5 = CGALi::hilbert_split (m4, m6, Cmp< z,  upz> (k));
+        Iterator m7 = CGALi::hilbert_split (m6, m8, Cmp< z, !upz> (k));
+
+        if (end - begin <= 8*8) return;
+
+        sort<z, upz, upx, upy> (m0, m1);
+        sort<y, upy, upz, upx> (m1, m2);
+        sort<y, upy, upz, upx> (m2, m3);
+        sort<x, upx,!upy,!upz> (m3, m4);
+        sort<x, upx,!upy,!upz> (m4, m5);
+        sort<y,!upy, upz,!upx> (m5, m6);
+        sort<y,!upy, upz,!upx> (m6, m7);
+        sort<z,!upz,!upx, upy> (m7, m8);
+    }
+
+    void operator() (Iterator begin, Iterator end) const
+    {
+        sort <0, false, false, false> (begin, end);
+    }
+};
+
+template <class RandomAccessIterator, class Kernel>
+void hilbert_sort_3 (RandomAccessIterator begin, RandomAccessIterator end,
+                     Kernel k)
+{
+    (Hilbert_sort_3<Kernel,RandomAccessIterator> (k)) (begin, end);
+}
+
+template <class RandomAccessIterator>
+void hilbert_sort_3 (RandomAccessIterator begin, RandomAccessIterator end)
+{
+    typedef std::iterator_traits<RandomAccessIterator> ITraits;
+    typedef typename ITraits::value_type               value_type;
+    typedef CGAL::Kernel_traits<value_type>            KTraits;
+    typedef typename KTraits::Kernel                   Kernel;
+
+    return hilbert_sort_3 (begin, end, Kernel());
+}
+
+CGAL_END_NAMESPACE
+
+#endif//CGAL_HILBERT_SORT_3_H
diff --git a/Spatial_sorting/include/CGAL/Hilbert_sort_base.h b/Spatial_sorting/include/CGAL/Hilbert_sort_base.h
new file mode 100644
index 00000000000..4043e71e50f
--- /dev/null
+++ b/Spatial_sorting/include/CGAL/Hilbert_sort_base.h
@@ -0,0 +1,27 @@
+#ifndef CGAL_HILBERT_SORT_BASE_H
+#define CGAL_HILBERT_SORT_BASE_H
+
+#include <CGAL/basic.h>
+#include <algorithm>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+    template <class RandomAccessIterator, class Cmp>
+    RandomAccessIterator
+    hilbert_split (RandomAccessIterator begin, RandomAccessIterator end,
+                   Cmp cmp = Cmp ())
+    {
+        if (begin >= end) return begin;
+
+        RandomAccessIterator middle = begin + (end - begin) / 2;
+        std::nth_element (begin, middle, end, cmp);
+        return middle;
+    }
+};
+
+CGAL_END_NAMESPACE
+
+
+#endif//CGAL_HILBERT_SORT_BASE_H
diff --git a/Spatial_sorting/maintainer b/Spatial_sorting/maintainer
new file mode 100644
index 00000000000..1f91c0fd21c
--- /dev/null
+++ b/Spatial_sorting/maintainer
@@ -0,0 +1 @@
+Christophe Delage <christophe.delage@sophia.inria.fr>
diff --git a/Spatial_sorting/test/Spatial_sorting/test_brio.C b/Spatial_sorting/test/Spatial_sorting/test_brio.C
new file mode 100644
index 00000000000..ff1fcde0a21
--- /dev/null
+++ b/Spatial_sorting/test/Spatial_sorting/test_brio.C
@@ -0,0 +1,89 @@
+#include <CGAL/basic.h>
+
+#include <CGAL/BRIO_sort.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/point_generators_2.h>
+#include <CGAL/point_generators_3.h>
+
+#include <iostream>
+#include <algorithm>
+#include <vector>
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_2 Point_2;
+typedef K::Point_3 Point_3;
+
+typedef CGAL::Creator_uniform_2<double,Point_2>  Creator_2;
+typedef CGAL::Creator_uniform_3<double,Point_3>  Creator_3;
+
+int main ()
+{
+    const int nb_points_2 = 50000, nb_points_3 = 50000;
+    CGAL::Random random (42);
+
+    std::cout << "Testing BRIO + Hilbert sort." << std::endl;
+
+    {
+        std::cout << "Testing 2D: Generating points... " << std::flush;
+
+        std::vector<Point_2> v;
+        v.reserve (nb_points_2);
+
+        CGAL::Random_points_in_square_2<Point_2,Creator_2> gen (1.0, random);
+
+        for (int i = 0; i < nb_points_2; ++i)
+            v.push_back (*gen++);
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Sorting points...    " << std::flush;
+
+        CGAL::brio_sort_2 (v.begin(), v.end());
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Checking...          " << std::flush;
+
+        std::vector<Point_2> v2 (v);
+        std::sort (v.begin(),  v.end(),  K().less_xy_2_object());
+        std::sort (v2.begin(), v2.end(), K().less_xy_2_object());
+        assert (v == v2);
+
+        std::cout << "OK." << std::endl;
+    }
+
+    {
+        std::cout << "Testing 3D: Generating points... " << std::flush;
+
+        std::vector<Point_3> v;
+        v.reserve (nb_points_3);
+
+        CGAL::Random_points_in_cube_3<Point_3,Creator_3> gen (1.0, random);
+
+        for (int i = 0; i < nb_points_3; ++i)
+            v.push_back (*gen++);
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Sorting points...    " << std::flush;
+
+        CGAL::brio_sort_3 (v.begin(), v.end());
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Checking...          " << std::flush;
+
+        std::vector<Point_3> v2 (v);
+        std::sort (v.begin(),  v.end(),  K().less_xyz_3_object());
+        std::sort (v2.begin(), v2.end(), K().less_xyz_3_object());
+        assert (v == v2);
+
+        std::cout << "OK." << std::endl;
+    }
+
+    return 0;
+}
diff --git a/Spatial_sorting/test/Spatial_sorting/test_hilbert.C b/Spatial_sorting/test/Spatial_sorting/test_hilbert.C
new file mode 100644
index 00000000000..cb04e9f515d
--- /dev/null
+++ b/Spatial_sorting/test/Spatial_sorting/test_hilbert.C
@@ -0,0 +1,90 @@
+#include <CGAL/basic.h>
+
+#include <CGAL/Hilbert_sort_2.h>
+#include <CGAL/Hilbert_sort_3.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/point_generators_2.h>
+#include <CGAL/point_generators_3.h>
+
+#include <iostream>
+#include <algorithm>
+#include <vector>
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_2 Point_2;
+typedef K::Point_3 Point_3;
+
+typedef CGAL::Creator_uniform_2<double,Point_2>  Creator_2;
+typedef CGAL::Creator_uniform_3<double,Point_3>  Creator_3;
+
+int main ()
+{
+    const int nb_points_2 = 50000, nb_points_3 = 50000;
+    CGAL::Random random (42);
+
+    std::cout << "Testing Hilbert sort." << std::endl;
+
+    {
+        std::cout << "Testing 2D: Generating points... " << std::flush;
+
+        std::vector<Point_2> v;
+        v.reserve (nb_points_2);
+
+        CGAL::Random_points_in_square_2<Point_2,Creator_2> gen (1.0, random);
+
+        for (int i = 0; i < nb_points_2; ++i)
+            v.push_back (*gen++);
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Sorting points...    " << std::flush;
+
+        CGAL::hilbert_sort_2 (v.begin(), v.end());
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Checking...          " << std::flush;
+
+        std::vector<Point_2> v2 (v);
+        std::sort (v.begin(),  v.end(),  K().less_xy_2_object());
+        std::sort (v2.begin(), v2.end(), K().less_xy_2_object());
+        assert (v == v2);
+
+        std::cout << "OK." << std::endl;
+    }
+
+    {
+        std::cout << "Testing 3D: Generating points... " << std::flush;
+
+        std::vector<Point_3> v;
+        v.reserve (nb_points_3);
+
+        CGAL::Random_points_in_cube_3<Point_3,Creator_3> gen (1.0, random);
+
+        for (int i = 0; i < nb_points_3; ++i)
+            v.push_back (*gen++);
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Sorting points...    " << std::flush;
+
+        CGAL::hilbert_sort_3 (v.begin(), v.end());
+
+        std::cout << "done." << std::endl;
+
+        std::cout << "            Checking...          " << std::flush;
+
+        std::vector<Point_3> v2 (v);
+        std::sort (v.begin(),  v.end(),  K().less_xyz_3_object());
+        std::sort (v2.begin(), v2.end(), K().less_xyz_3_object());
+        assert (v == v2);
+
+        std::cout << "OK." << std::endl;
+    }
+
+    return 0;
+}