Edit on GitHub

Usage

All of the Hilbert algorithms have the same interface.

Decide dimensions of the spatial axes

All of the algorithms, except Simple2D, need to know ahead of time the extent of the coordinate system. If the Hilbert curve will be in a 32 x 32 space, then the dimension is 2, and it needs log2(32)=5 bits of resolution in both directions. For GlobalGray, that's b=5, n=2. If the sizes aren't powers of two or are uneven, then set the bits to cover the largest side, so (12 x 12 x 12 x 800) would be b=10, n=4 because $800 < 2^{10}$. There is one algorithm that deals better with uneven sides. The Compact algorithm can pack bits together so that the resulting Hilbert index takes less storage space. This is usually used for database storage. It could take (12 x 12 x 12 x 800) as Compact([4, 4, 4, 10]) which results in a 24-bit integer that can be stored in a UInt32.

Create an algorithm

For most Hilbert curve algorithms, you have to say, beforehand, how large the multidimensional coordinates are, in powers of two. For instance, a three-dimensional grid can have values from 1 to 16 in each dimension, so n = 3 and b = 4 because 16 = 2^b.

using BijectiveHilbert
dimensions = 3
bits = 4
simple = Simple2D(Int)
gg = GlobalGray(UInt, bits, dimensions)
sg = SpaceGray(UInt, bits, dimensions)
compact = Compact(UInt, fill(bits, dimensions))

The first argument is a datatype for the Hilbert index. It should be large enough to hold all of the bits from the n-dimensional axes. If you don't specify one, it will use the smallest unsigned integer that can hold them.

Note that the Compact algorithm can have different sizes for each dimension, but they are all powers of two. It produces a Hilbert index that uses only as many bits as necessary.

Encode and decode

You can encode from n-dimensions to the Hilbert index, or you can decode from a Hilbert index to n-dimensions.

for algorithm in [simple, gg, sg, compact]
    for k in 1:(1<<bits)
        for j in 1:(1<<bits)
            for i in 1:(1<<bits)
                X = [i, j, k]
                h = encode_hilbert(algorithm, X)
                X .= 0
                decode_hilbert!(algorithm, X, h)
            end
        end
    end
end

Encode and decode with a zero-based value

The underlying algorithms use a zero-based axis and a zero-based Hilbert index. These are available, too.

for algorithm in [simple, gg, sg, compact]
    for k in 0:(1<<bits - 1)
        for j in 0:(1<<bits - 1)
            for i in 0:(1<<bits - 1)
                X = [i, j, k]
                h = encode_hilbert_zero(algorithm, X)
                X .= 0
                decode_hilbert_zero!(algorithm, X, h)
            end
        end
    end
end

Index

BijectiveHilbert.Bijective.Simple2D — Type

Simple2D(::Type{T})

The type is a data type to hold the Hilbert index. It should have enough bits to hold all bits of integer axes to encode.

The variant algorithm used differs from the usual Hilbert code because it doesn't need to know the size of the whole grid before computing the code. It looks like a slightly-rotated version of the Hilbert curve, but it has the benefit that it is 1-1 between (x, y) and z, so you can translate back and forth.

If the size of the axis values or the size of the Hilbert index are too large to be stored in the data types of the axis vector or index, then you'll see an error from a failure of trunc, which tries to copy values into the smaller datatype. Solve this by using a larger datatype, so UInt64 instead of UInt32.

It comes from a paper: N. Chen, N. Wang, B. Shi, A new algorithm for encoding and decoding the Hilbert order. Software—Practice and Experience 2007; 37(8): 897–908.

source
BijectiveHilbert.Compact — Type
Compact(ms::Vector{Int})
Compact(::Type{T}, ms::Vector{Int})

This algorithm is n-dimensional and permits dimensions to use different numbers of bits, specified in the ms vector. The type T is an optional data type for the Hilbert index. It should be greater than or equal to the sum of the bits.

This algorithm comes from three sources:

  • A technical report, "Compact Hilbert Indices" by Chris Hamilton. Technical Report CS-2006-07. 6059 University Ave., Halifax, Nova Scotia, B3H 1W5, Canada. This report is informative but has many errors.

  • A paper by Hamilton and Rau-Chaplin, "Compact Hilbert Indices for Multi-Dimensional Data," 2007. Nice paper. Also wrong.

  • The libhilbert source code is a copy of Hamilton's work and has many corrections. This, ultimately, lead to the working code.

source
BijectiveHilbert.FaceContinuous — Type
FaceContinuous(b::Int, n::Int)
FaceContinuous(::Type{T}, b::Int, n::Int)

For n dimensions, use b bits of precision in this Hilbert curve. If you specify a type T, this will be used as the type of the Hilbert encoding. If not, the smallest unsigned integer that can hold n*b bits will be used for the Hilbert index data type.

This is the Butz algorithm, as presented by Lawder. Haverkort2017 says it is face continuous. The code is in lawder.c. The original paper had an error, and Lawder put a correction on his website. http://www.dcs.bbk.ac.uk/~jkl/publications.html

source
BijectiveHilbert.GlobalGray — Type
GlobalGray(b, n)
GlobalGray(T, b, n)

T is a data type for the Hilbert index. It can be signed or unsigned, as long as it has at least n * b bits. n is the number of dimensions, and b is the bits per dimension, so each axis value should be between 0 and $2^b - 1$, inclusive, for the zero-based interface. They should be between 1 and $2^b$, inclusive, for the one-based interface.

The GlobalGray algorithm is an n-dimensional Hilbert curve with a simplified implementation. It follows an article, "Programming the Hilbert Curve," by John Skilling, 707 (2004), http://dx.doi.org/10.1063/1.1751381. I call it "Global Gray" because the insight of the article is that a single, global Gray code can be applied to all np bits of a Hilbert length.

source
BijectiveHilbert.SpaceGray — Type
SpaceGray(b, n)
SpaceGray(::Type{T}, b, n)

This is an n-dimensional Hilbert curve where all n dimensions must have b bits in size. It was described in the same paper and examples as the Compact algorithm.

source
BijectiveHilbert.decode_hilbert_zero! — Method
decode_hilbert_zero!(ha::HilbertAlgorithm{T}}, X::Vector{A}, h::T)

Given a Hilbert index, h, computes an n-dimensional coordinate X. The type of the Hilbert index is large enought to contain the bits of all dimensions of the axis vector, X.

source
BijectiveHilbert.encode_hilbert — Method
encode_hilbert(ha::HilbertAlgorithm{T}, X::Vector{A})

A 1-based Hilbert encoding. Both the Hilbert index and the axes start counting at 1 instead of 0.

source
BijectiveHilbert.encode_hilbert_zero — Method
encode_hilbert_zero(ha::HilbertAlgorithm{T}, X::Vector{A})

Takes an n-dimensional vector X and returns a single integer of type T which orders X to improve spatial locality. The input X has multiple axes and the output is called a Hilbert index. This version is zero-based, so each axis counts from 0, and the smallest Hilbert index is 0.

source