The Hilbert curve is a fractal space-filling curve
that is rather pretty to look at.
A Hilbert curve of order n
traces a single path over a square of side 2^n
units, as you can see in the images from MathWorld above (with curves of order 2 through 6).
If you pick a point p
from within the square covered by a curve, and trace out the curve until you reach that point, the number of intermediate points that you passed through along your way to p
is called the Hilbert distance. The notion of distance along a space-filling curve is useful for a few reasons.
- You can use it as a hash from a 2D coordinate to an integer.
- Some curves preserve locality; for example, points that are nearby in 2D space are likely to have similar Hilbert distances.
- It’s an isomorphism: you can turn a distance back into a 2D coordinate, in some cases quite easily.
- Locality-preserving hashes can be useful for clustering nodes of an R-tree, which is the generalisation of a B-tree to spatial data.
There’s a nice algorithm for computing the Hilbert distance of a 2D point that’s linear in the order of the curve (or logarithmic in the size of the side of the square you’re trying to cover). To put this into useful terms, a “10 megapixel” camera typically has a resolution of 3766×2800 pixels. To cover every pixel with a Hilbert curve, you’d need a curve with 4096 points on a side, which is an order-12 curve. You can compute the Hilbert distance of any point in your 10-megapixel photo with 12 steps of the distance algorithm.
Here’s an implementation of this Hilbert distance algorithm in Haskell.
hilbertDistance :: (Bits a, Ord a) => Int -> (a,a) -> a
hilbertDistance d (x,y)
| x < 0 || x >= 1 `shiftL` d = error "x bounds"
| y < 0 || y >= 1 `shiftL` d = error "y bounds"
| otherwise = dist (1 `shiftL` (d - 1)) (1 `shiftL` ((d - 1) * 2)) 0 x y
where dist 0 _ result _ _ = result
dist side area result x y =
case (compare x side, compare y side) of
(LT, LT) -> step result y x
(LT, _) -> step (result + area) x (y - side)
(_, LT) -> step (result + area * 3) (side - y - 1)
(side * 2 - x - 1)
(_, _) -> step (result + area * 2) (x - side) (y - side)
where step = dist (side `shiftR` 1) (area `shiftR` 2)
As you can see, each step requires a few comparisons, a few shifts, and a handful of arithmetic operations, so it’s fairly cheap, and the loop unrolls and constant-folds well.
There’s some controversy in the literature over whether a Hilbert curve gives better spatial locality than a Lebesgue (also known as a Z-order, or Morton) curve. It’s certainly the case that computing the Lebesgue distance is much cheaper than the Hilbert distance, as it’s possible to finagle an answer using a fixed, unconditional sequence of bit-twiddling operations
to interleave the bits of the x
coordinate with those of the y
lebesgueDistance :: (Int, Int) -> Int
lebesgueDistance (x,y) =
let ms = [(0x00FF00FF, 8), (0x0F0F0F0F, 4),
(0x33333333, 2), (0x55555555, 1)]
aaargh x (mask,shift) = (x .|. (x `shiftL` shift)) .&. mask
x' = foldl aaargh x ms
y' = foldl aaargh y ms
in x .|. (y `shiftL` 1)
I haven’t myself checked to see which curve better preserves locality. It’s certainly easier to see how to extend the Lebesgue bit-fiddling hack to a higher dimension, while the given Hilbert code is a bit of a head wrecker.