/ 
Math-PlanePath-13
River stage one • 3 direct dependents • 4 total dependents
/ Math::PlanePath::HilbertCurve

NAME

Math::PlanePath::HilbertCurve -- self-similar quadrant traversal

SYNOPSIS

use Math::PlanePath::HilbertCurve;
my $path = Math::PlanePath::HilbertCurve->new;
my ($x, $y) = $path->n_to_xy (123);

DESCRIPTION

This path by David Hilbert traverses a quadrant of the plane one step at a time in a self-similar pattern,

       ...
        |
y=7    63--62  49--48--47  44--43--42
            |   |       |   |       |
y=6    60--61  50--51  46--45  40--41
        |           |           |
y=5    59  56--55  52  33--34  39--38
        |   |   |   |   |   |       |
y=4    58--57  54--53  32  35--36--37
                        |
y=3     5---6   9--10  31  28--27--26
        |   |   |   |   |   |       |
y=2     4   7---8  11  30--29  24--25
        |           |           |
y=1     3---2  13--12  17--18  23--22
            |   |       |   |       |
y=0     0---1  14--15--16  19--20--21

      x=0   1   2   3   4   5   6   7

The start is a sideways U shape per 0,1,2,3, and then four of those are put together in an upside-down U. The orientation of the sub parts are chosen so the starts and ends are adjacent, so 3 next to 4, 7 next to 8, and 11 next to 12.

5,6___9,10
4,7   8,11
 |     |
3,2   13,12__
0,1   14,15

The process repeats, doubling in size each time and alternately sideways or upside-down U at the top level and invert or transponses as necessary in the sub-parts.

The pattern can be drawn with the first step 0->1 up instead of to the right. Right is used here since that's what most of the other PlanePaths do. Swap X and Y for upwards first instead.

Within a power-of-2 square 2x2, 4x4, 8x8, 16x16 etc 2^(2^k), all the N values 0 to 2^(2*(2^k))-1 are within the square. The alternate top left or bottom right corner 3, 15, 63, 255 etc of each is the 2^(2*(2^k))-1 maximum.

OEIS

The Hilbert Curve path is in Sloane's OEIS as sequences A163355, A163357, A163359 and A163361. They're based on numbering X,Y positions diagonally in the style of Math::PlanePath::Diagonals, but starting from N=0. The four sequences are whether the first Curve move is in the X or Y direction, and then whether the diagonals are numbered from the X axis up to the Y or the other way around.

The sequences are permutations of the integers, and A163356, A163358, A163360 and A163362 are the corresponding inverses.

Algorithms

Converting an N to X,Y coordinates is reasonably straightforward. The top two bits of N is a configuration

3--2                    1--2
   |    or transpose    |  |
0--1                    0  3

according to whether it's an odd or even bit-pair position. Within the "3" sub-parts there's also inverted forms

1--0        3  0
|           |  |
2--3        2--1

Working N from high to low a state variable can record whether there's a transpose, an invert, or both (four states altogether). A bit pair 0,1,2,3 from N then gives a bit each of X,Y according to the configuration, and a new state which is the orientation of the sub-part.

Gosper's HAKMEM item 115 has this with tables for the state and X,Y bits,

http://www.inwap.com/pdp10/hbaker/hakmem/topology.html#item115

And C++ code based on that in Jorg Arndt's book,

http://www.jjj.de/fxt/#fxtbook (section 1.31.1)

It also works to process N from low to high, at each stage applying a transpose (swap X,Y) and/or invert (bitwise negate) to the low X,Y bits generated so far. This approach saves locating the top bits of N, but if using bignums then the bitwise inverts will be much more work.

The reverse X,Y to N can follow the table approach from high to low taking one bit from X and Y each time. The state-based table of N-pair -> X-bit,Y-bit is reversible and the new state is based on the N-pair so obtained (or based on the X,Y bits if that mapping is combined in).

The current code is a mixture of low to high for n_to_xy but the table high to low for the reverse xy_to_n.

Each step between successive N values is by 1 up, down, left or right. The next direction can be calculated from the N position based on some base-4 parity of N and -N. C++ code in Jorg Arndt's fxtbook per above.

FUNCTIONS

$path = Math::PlanePath::HilbertCurve->new ()

Create and return a new path object.

($x,$y) = $path->n_to_xy ($n)

Return the X,Y coordinates of point number $n on the path. Points begin at 0 and if $n < 0 then the return is an empty list.

Fractional positions give an X,Y position along a straight line between the integer positions. Integer positions are always just 1 apart either horizontally or vertically, so the effect is that the fraction part appears either added to or subtracted from X or Y.

$n = $path->xy_to_n ($x,$y)

Return an integer point number for coordinates $x,$y. Each integer N is considered the centre of a unit square an $x,$y within that square returns N.

SEE ALSO

Math::PlanePath, Math::PlanePath::ZOrderCurve

HOME PAGE

http://user42.tuxfamily.org/math-planepath/index.html

LICENSE

Math-PlanePath is Copyright 2010 Kevin Ryde

Math-PlanePath is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version.

Math-PlanePath is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.