diff --git a/include/geode/geometry/points_sort.hpp b/include/geode/geometry/points_sort.hpp
index 074a6350c..c1eeaa574 100644
--- a/include/geode/geometry/points_sort.hpp
+++ b/include/geode/geometry/points_sort.hpp
@@ -42,4 +42,8 @@ namespace geode
     template < index_t dimension >
     [[nodiscard]] std::vector< index_t > morton_mapping(
         absl::Span< const Point< dimension > > points );
+
+    template < index_t dimension >
+    [[nodiscard]] std::vector< index_t > hilbert_mapping(
+        absl::Span< const Point< dimension > > points );
 } // namespace geode
diff --git a/src/geode/geometry/points_sort.cpp b/src/geode/geometry/points_sort.cpp
index b6ae4829c..d2c099112 100644
--- a/src/geode/geometry/points_sort.cpp
+++ b/src/geode/geometry/points_sort.cpp
@@ -57,6 +57,28 @@ namespace
     };
     ALIAS_1D_AND_2D_AND_3D( Morton_cmp );
 
+    template < geode::index_t dimension >
+    class Hilbert_cmp
+    {
+    public:
+        Hilbert_cmp( absl::Span< const geode::Point< dimension > > points,
+            geode::local_index_t coord )
+            : points_( points ), coord_( coord )
+        {
+        }
+
+        bool operator()( geode::index_t box1, geode::index_t box2 ) const
+        {
+            return points_[box1].value( coord_ )
+                   > points_[box2].value( coord_ );
+        }
+
+    private:
+        absl::Span< const geode::Point< dimension > > points_;
+        geode::local_index_t coord_;
+    };
+    ALIAS_1D_AND_2D_AND_3D( Hilbert_cmp );
+
     /**
      * \brief Splits a sequence into two ordered halves.
      * \details The algorithm shuffles the sequence and
@@ -88,7 +110,8 @@ namespace
      *   In CGAL User and Reference Manual. CGAL Editorial Board,
      *   3.9 edition, 2011
      */
-    template < geode::local_index_t COORDX >
+    template < geode::local_index_t COORDX,
+        template < geode::index_t > typename Comparator >
     void morton_mapping( absl::Span< const geode::Point3D > points,
         const itr& begin,
         const itr& end )
@@ -100,9 +123,9 @@ namespace
         constexpr geode::local_index_t COORDY = COORDX == 2 ? 0 : COORDX + 1;
         constexpr geode::local_index_t COORDZ = COORDY == 2 ? 0 : COORDY + 1;
 
-        const Morton_cmp3D compX{ points, COORDX };
-        const Morton_cmp3D compY{ points, COORDY };
-        const Morton_cmp3D compZ{ points, COORDZ };
+        const Comparator< 3 > compX{ points, COORDX };
+        const Comparator< 3 > compY{ points, COORDY };
+        const Comparator< 3 > compZ{ points, COORDZ };
 
         const auto m0 = begin;
         const auto m8 = end;
@@ -113,17 +136,18 @@ namespace
         const auto m6 = split_container( m4, m8, compY );
         const auto m5 = split_container( m4, m6, compZ );
         const auto m7 = split_container( m6, m8, compZ );
-        morton_mapping< COORDZ >( points, m0, m1 );
-        morton_mapping< COORDY >( points, m1, m2 );
-        morton_mapping< COORDY >( points, m2, m3 );
-        morton_mapping< COORDX >( points, m3, m4 );
-        morton_mapping< COORDX >( points, m4, m5 );
-        morton_mapping< COORDY >( points, m5, m6 );
-        morton_mapping< COORDY >( points, m6, m7 );
-        morton_mapping< COORDZ >( points, m7, m8 );
+        morton_mapping< COORDZ, Comparator >( points, m0, m1 );
+        morton_mapping< COORDY, Comparator >( points, m1, m2 );
+        morton_mapping< COORDY, Comparator >( points, m2, m3 );
+        morton_mapping< COORDX, Comparator >( points, m3, m4 );
+        morton_mapping< COORDX, Comparator >( points, m4, m5 );
+        morton_mapping< COORDY, Comparator >( points, m5, m6 );
+        morton_mapping< COORDY, Comparator >( points, m6, m7 );
+        morton_mapping< COORDZ, Comparator >( points, m7, m8 );
     }
 
-    template < geode::local_index_t COORDX >
+    template < geode::local_index_t COORDX,
+        template < geode::index_t > typename Comparator >
     void morton_mapping( absl::Span< const geode::Point2D > points,
         const itr& begin,
         const itr& end )
@@ -134,21 +158,22 @@ namespace
         }
         constexpr geode::local_index_t COORDY = COORDX == 1 ? 0 : COORDX + 1;
 
-        const Morton_cmp2D compX{ points, COORDX };
-        const Morton_cmp2D compY{ points, COORDY };
+        const Comparator< 2 > compX{ points, COORDX };
+        const Comparator< 2 > compY{ points, COORDY };
 
         const auto m0 = begin;
         const auto m4 = end;
         const auto m2 = split_container( m0, m4, compX );
         const auto m1 = split_container( m0, m2, compY );
         const auto m3 = split_container( m2, m4, compY );
-        morton_mapping< COORDY >( points, m0, m1 );
-        morton_mapping< COORDX >( points, m1, m2 );
-        morton_mapping< COORDX >( points, m2, m3 );
-        morton_mapping< COORDY >( points, m3, m4 );
+        morton_mapping< COORDY, Comparator >( points, m0, m1 );
+        morton_mapping< COORDX, Comparator >( points, m1, m2 );
+        morton_mapping< COORDX, Comparator >( points, m2, m3 );
+        morton_mapping< COORDY, Comparator >( points, m3, m4 );
     }
 
-    template < geode::local_index_t COORDX >
+    template < geode::local_index_t COORDX,
+        template < geode::index_t > typename Comparator >
     void morton_mapping( absl::Span< const geode::Point1D > points,
         const itr& begin,
         const itr& end )
@@ -157,13 +182,13 @@ namespace
         {
             return;
         }
-        const Morton_cmp1D compX{ points, COORDX };
+        const Comparator< 1 > compX{ points, COORDX };
 
         const auto m0 = begin;
         const auto m2 = end;
         const auto m1 = split_container( m0, m2, compX );
-        morton_mapping< COORDX >( points, m0, m1 );
-        morton_mapping< COORDX >( points, m1, m2 );
+        morton_mapping< COORDX, Comparator >( points, m0, m1 );
+        morton_mapping< COORDX, Comparator >( points, m1, m2 );
     }
 
     /*
@@ -214,7 +239,22 @@ namespace geode
             [&mapping]( index_t i ) {
                 mapping[i] = i;
             } );
-        ::morton_mapping< 0_uc >( points, mapping.begin(), mapping.end() );
+        ::morton_mapping< 0_uc, Morton_cmp >(
+            points, mapping.begin(), mapping.end() );
+        return mapping;
+    }
+
+    template < index_t dimension >
+    std::vector< index_t > hilbert_mapping(
+        absl::Span< const Point< dimension > > points )
+    {
+        std::vector< index_t > mapping( points.size() );
+        async::parallel_for( async::irange( size_t{ 0 }, mapping.size() ),
+            [&mapping]( index_t i ) {
+                mapping[i] = i;
+            } );
+        ::morton_mapping< 0_uc, Hilbert_cmp >(
+            points, mapping.begin(), mapping.end() );
         return mapping;
     }
 
@@ -231,4 +271,11 @@ namespace geode
         absl::Span< const Point< 2 > > );
     template std::vector< index_t > opengeode_geometry_api morton_mapping(
         absl::Span< const Point< 3 > > );
+
+    template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
+        absl::Span< const Point< 1 > > );
+    template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
+        absl::Span< const Point< 2 > > );
+    template std::vector< index_t > opengeode_geometry_api hilbert_mapping(
+        absl::Span< const Point< 3 > > );
 } // namespace geode