Hilbert Curve

I’ve been reading a bit about the spacefilling curves for my wavelet image compression (take a look here and here).
There is a very nice way to convert from the Hilbert derived key to (multi-dimensional coordinates) described by John J. Bartholdi, III and Paul Goldsman in “Vertex-Labeling Algorithms for the Hilbert Spacefilling Curve”.
They describe a recursive procedure, but in the particular case mentioned above, this can be easily unrolled. It also works very well with fixed point arithmetic. The following source code can be further optimized by storing each point’s x and y coordinate in a single unsigned int, as everything except for the final averaging is easily replicated across the vector by applying the operation to the combined integer.

/** Converts a Hilbert Curve derived key into the coordinates of a point.
* This implementation is based on the paper "Vertex-Labeling Algorithms for the
* Hilbert Spacefilling Curve" by John J. Bartholdi, III and Paul Goldsman. They
* describe a recursive procedure which has been unrolled here. Also, due to
* fixed point arithmetic, the coordinates of the points grow with each
* iteration so that sufficient accuracy is maintained.
* @param[in] key Hilbert derived key to be converted. In essence, we interpret
*   every two bits of the key as a certain spatial subdivision of the
*   unit-square, as described in the paper.
* @param[in] numbits number of bits in the key (which has to be even). Only
*   valid for even values smaller than or equal to 32 (as otherwise key does
*   not have enough bits).
* @param[out] xx fixed point x-coordinate with a precision of numbits / 2.
* @param[out] yy fixed point y-coordinate with a precision of numbits / 2.
**/
void convert_hilbert_key(unsigned int key, int numbits, unsigned int* xx, unsigned int* yy)
{
unsigned int abcd[4][2] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; /* unit square */

while (numbits > 1) /* should be > 0, but this is safer for (invalid) odd numbers */
{
unsigned int tmp[4][2]; /* save abcd here */
unsigned char subcell;

memcpy(tmp, abcd, sizeof tmp); /* save our 4 points with the new ones */

numbits -= 2; /* each subdivision decision takes 2 bits */
subcell = (key >> numbits) & 3; /* namely these two (taken from the top) */
key &= (1 < < numbits) - 1; /* update key by removing the bits we used */

/** Add the two points with indices u and v (in tmp) and store the result at
* index dst in abcd (for both x(0) and y(1)).
**/
#define ADD(dst, u, v) (abcd[(dst)][0] = tmp[(u)][0] + tmp[(v)][0],
abcd[(dst)][1] = tmp[(u)][1] + tmp[(v)][1])

switch (subcell)
{ /* divide into subcells */
case 0:
/* h(key, numbits, a << 1, a + d, a + c, a + b); */
ADD(0, 0, 0); ADD(1, 0, 3); ADD(2, 0, 2); ADD(3, 0, 1);
break;
case 1:
/* h(key, numbits, b + a, b << 1, b + c, b + d); */
ADD(0, 1, 0); ADD(1, 1, 1); ADD(2, 1, 2); ADD(3, 1, 3);
break;
case 2:
/* h(key, numbits, c + a, c + b, c << 1, c + d); */
ADD(0, 2, 0); ADD(1, 2, 1); ADD(2, 2, 2); ADD(3, 2, 3);
break;
case 3:
/* h(key, numbits, d + c, d + b, d + a, d << 1); */
ADD(0, 3, 2); ADD(1, 3, 1); ADD(2, 3, 0); ADD(3, 3, 3);
break;
}

#undef ADD
}
/* final result is the midpoint of the cell, i.e. (a + b + c + d) / 4 */
*xx = (abcd[0][0] + abcd[1][0] + abcd[2][0] + abcd[3][0] + 1) >> 2;
*yy = (abcd[0][1] + abcd[1][1] + abcd[2][1] + abcd[3][1] + 1) >> 2;
/*	printf("x: %u y: %un", *xx, *yy);*/
}

Leave a Reply

Your email address will not be published.