1Math::PlanePath::HilberUtsCeurrvCeo(n3t)ributed Perl DocMuamtehn:t:aPtliaonnePath::HilbertCurve(3)
2
3
4

NAME

6       Math::PlanePath::HilbertCurve -- 2x2 self-similar quadrant traversal
7

SYNOPSIS

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

DESCRIPTION

14       This path is an integer version of the curve described by Hilbert in
15       1891 for filling a unit square.  It traverses a quadrant of the plane
16       one step at a time in a self-similar 2x2 pattern,
17
18           David Hilbert, "Ueber die stetige Abbildung einer Linie auf ein
19           Flächenstück", Mathematische Annalen, volume 38, number 3, 1891,
20           pages 459-460, DOI 10.1007/BF01199431.
21
22                  ...
23               |   |
24             7 |  63--62  49--48--47  44--43--42
25               |       |   |       |   |       |
26             6 |  60--61  50--51  46--45  40--41
27               |   |           |           |
28             5 |  59  56--55  52  33--34  39--38
29               |   |   |   |   |   |   |       |
30             4 |  58--57  54--53  32  35--36--37
31               |                   |
32             3 |   5---6   9--10  31  28--27--26
33               |   |   |   |   |   |   |       |
34             2 |   4   7---8  11  30--29  24--25
35               |   |           |           |
36             1 |   3---2  13--12  17--18  23--22
37               |       |   |       |   |       |
38           Y=0 |   0---1  14--15--16  19--20--21
39               +----------------------------------
40                 X=0   1   2   3   4   5   6   7
41
42       The start is a sideways U shape N=0 to N=3, then four of those are put
43       together in an upside-down U as
44
45           5,6    9,10
46           4,7--- 8,11
47             |      |
48           3,2   13,12
49           0,1   14,15--
50
51       The orientation of the sub parts ensure the starts and ends are
52       adjacent, so 3 next to 4, 7 next to 8, and 11 next to 12.
53
54       The process repeats, doubling in size each time and alternately
55       sideways or upside-down U with invert and/or transpose as necessary in
56       the sub-parts.
57
58       The pattern is sometimes drawn with the first step 0->1 upwards instead
59       of to the right.  First step right is used here for consistency with
60       other PlanePaths.  Swap X and Y for upwards first instead.
61
62       See examples/hilbert-path.pl for a sample program printing the path
63       pattern in ascii.
64
65   Level Ranges
66       Within a power-of-2 square 2x2, 4x4, 8x8, 16x16 etc (2^k)x(2^k) at the
67       origin, all the N values 0 to 2^(2*k)-1 are within the square.  The
68       maximum 3, 15, 63, 255 etc 2^(2*k)-1 is alternately at the top left or
69       bottom right corner.
70
71       Because each step is by 1, the distance along the curve between two X,Y
72       points is the difference in their N values (as from xy_to_n()).
73
74       On the X=Y diagonal N=0,2,8,10,32,etc is the integers using only digits
75       0 and 2 in base 4, or equivalently have even-numbered bits 0, like
76       x0y0...z0.
77
78   Locality
79       The Hilbert curve is fairly well localized in the sense that a small
80       rectangle (or other shape) is usually a small range of N.  This
81       property is used in some database systems to store X,Y coordinates
82       using the resulting Hilbert curve N as an index.  A search through a
83       2-D region is then usually a fairly modest linear search through N
84       values.  rect_to_n_range() gives exact N range for a rectangle, or see
85       notes "Rectangle to N Range" below for calculating on any shape.
86
87       The N range can be large when crossing sub-parts.  In the sample above
88       it can be seen for instance adjacent points X=0,Y=3 and X=0,Y=4 have
89       rather widely spaced N values 5 and 58.
90
91       Fractional X,Y values can be indexed by extending the N calculation
92       down into X,Y binary fractions.  The code here doesn't do that, but
93       could be pressed into service by moving the binary point in X and Y an
94       even number of places, the same in each.  (An odd number of bits would
95       require swapping X,Y to compensate for the alternating transpose in
96       part 0.)  The resulting integer N is then divided down by a
97       corresponding multiple-of-4 binary places.
98

FUNCTIONS

100       See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
101       classes.
102
103       "$path = Math::PlanePath::HilbertCurve->new ()"
104           Create and return a new path object.
105
106       "($x,$y) = $path->n_to_xy ($n)"
107           Return the X,Y coordinates of point number $n on the path.  Points
108           begin at 0 and if "$n < 0" then the return is an empty list.
109
110           Fractional positions give an X,Y position along a straight line
111           between the integer positions.  Integer positions are always just 1
112           apart either horizontally or vertically, so the effect is that the
113           fraction part is an offset along either $x or $y.
114
115       "$n = $path->xy_to_n ($x,$y)"
116           Return an integer point number for coordinates "$x,$y".  Each
117           integer N is considered the centre of a unit square and an "$x,$y"
118           within that square returns N.
119
120       "($n_lo, $n_hi) = $path->rect_to_n_range ($x1,$y1, $x2,$y2)"
121           The returned range is exact, meaning $n_lo and $n_hi are the
122           smallest and biggest in the rectangle.
123
124   Level Methods
125       "($n_lo, $n_hi) = $path->level_to_n_range($level)"
126           Return "(0, 4**$level - 1)".
127

FORMULAS

129   N to X,Y
130       Converting N to X,Y coordinates is reasonably straightforward.  The top
131       two bits of N is a configuration
132
133           3--2                    1--2
134              |    or transpose    |  |
135           0--1                    0  3
136
137       according to whether it's an odd or even bit-pair position.  Then
138       within each of the "3" sub-parts there's also inverted forms
139
140           1--0        3  0
141           |           |  |
142           2--3        2--1
143
144       Working N from high to low with a state variable can record whether
145       there's a transpose, an invert, or both, being four states altogether.
146       A bit pair 0,1,2,3 from N then gives a bit each of X,Y according to the
147       configuration and a new state which is the orientation of that sub-
148       part.  Bill Gosper's HAKMEM item 115 has this with either bit
149       operations or a table for the state and X,Y bits,
150
151           <https://dspace.mit.edu/handle/1721.1/6086>,
152           <http://www.inwap.com/pdp10/hbaker/hakmem/topology.html#item115>
153
154       And C++ code based on that in Jorg Arndt's book,
155
156           <http://www.jjj.de/fxt/#fxtbook> (section 1.31.1)
157
158       It also works to process N from low to high, at each stage applying any
159       transpose (swap X,Y) and/or invert (bitwise NOT) to the low X,Y bits
160       generated so far.  This works because there's no "reverse" sections, or
161       equivalently the curve is the same forward and reverse.  Low to high
162       saves locating the top bits of N, but if using bignums then the bitwise
163       inverts of the full X,Y values will be much more work.
164
165   X,Y to N
166       X,Y to N can follow the table approach from high to low taking one bit
167       from X and Y each time.  The state table of N-pair -> X-bit,Y-bit is
168       reversible, and a new state is based on the N-pair thus obtained (or
169       could be based on the X,Y bits if that mapping is combined into the
170       state transition table).
171
172   Rectangle to N Range
173       An easy over-estimate of the maximum N in a region can be had by
174       finding the next bigger (2^k)x(2^k) square enclosing the region.  This
175       means the biggest X or Y rounded up to the next power of 2, so
176
177           find lowest k with 2^k > max(X,Y)
178           N_max = 2^(2k) - 1
179
180       Or equivalently rounding down to the next lower power of 2,
181
182           find highest k with 2^k <= max(X,Y)
183           N_max = 2^(2*(k+1)) - 1
184
185       An exact N range can be found by following the high to low N to X,Y
186       procedure above.  Start at the 2^(2k) bit pair position in an N bigger
187       than the desired region and choose 2 bits for N to give a bit each of X
188       and Y.  The X,Y bits are based on the state table as above and the bits
189       chosen for N are those for which the resulting X,Y sub-square overlaps
190       some of the target region.  The smallest N similarly, choosing the
191       smallest bit pair for N which overlaps.
192
193       The biggest and smallest N digit for a sub-part can be found with a
194       lookup table.  The X range might cover one or both sub-parts, and the Y
195       range similarly, for a total 4 possible configurations.  Then a table
196       of state+coverage -> digit gives the minimum and maximum N bit-pair,
197       and state+digit gives a new state the same as X,Y to N.
198
199       Biggest and smallest N must be calculated with separate state and X,Y
200       values since they track down different N bits and thus different
201       states.  But they take the same number of steps from an enclosing level
202       down to level 0 and can thus be done in a single loop.
203
204       The N range for any shape can be found this way, not just a rectangle
205       like rect_to_n_range().  At each level the procedure only depends on
206       asking which of the four sub-parts overlaps some of the target area.
207
208   Direction
209       Each step between successive N values is always 1 up, down, left or
210       right.  The next direction can be calculated from N in the high-to-low
211       procedure above by watching for the lowest non-3 digit and noting the
212       direction from that digit towards digit+1.  That can be had from the
213       state+digit -> X,Y table looking up digit and digit+1, or alternatively
214       a further table encoding state+digit -> direction.
215
216       The reason for taking only the lowest non-3 digit is that in a 3 sub-
217       part the direction it goes is determined by the next higher level.  For
218       example at N=11 the direction is down for the inverted-U of the next
219       higher level N=0,4,8,12.
220
221       This non-3 (or non whatever highest digit) is a general procedure and
222       can be used on any state-based high-to-low procedure of self-similar
223       curves.  In the current code, it's used to apply a fractional part of N
224       in the correct direction but is not otherwise made directly available.
225
226       Because the Hilbert curve has no "reverse" sections it also works to
227       build a direction from low to high N digits.  1 and 2 digits make no
228       change to the orientation, 0 digit is a transpose, and a 3 digit is a
229       rotate and transpose, except that low 3s are transpose-only (no rotate)
230       for the same reason as taking the lowest non-3 above.
231
232       Jorg Arndt in the fxtbook above notes the direction can be obtained
233       just by counting 3s in n and -n (the twos-complement).  This numbers
234       segments starting n=1, unlike PlanePath here starting N=0, so it
235       becomes
236
237           N+1 count 3s  / 0 mod 2   S or E
238                         \ 1 mod 2   N or W
239
240           -(N+1) count 3s  / 0 mod 2   N or E
241                            \ 1 mod 2   S or W
242
243       For the twos-complement negation, an even number of base-4 digits of N
244       must be taken.  Because -(N+1) = ~N, ie. a ones-complement, the second
245       part is also
246
247           N count 0s          / 0 mod 2   N or E
248           in even num digits  \ 1 mod 2   S or W
249
250       Putting the two together then
251
252           N count 0s   N+1 count 3s    direction (0=E,1=N,etc)
253           in base 4    in base 4
254
255             0 mod 2      0 mod 2          0
256             1 mod 2      0 mod 2          3
257             0 mod 2      1 mod 2          1
258             1 mod 2      1 mod 2          2
259
260   Segments in Direction
261       The number of segments in each direction is calculated in
262
263           Sergey Kitaev, Toufik Mansour and Patrice Séébold, "Generating the
264           Peano Curve and Counting Occurrences of Some Patterns", Journal of
265           Automata, Languages and Combinatorics, volume 9, number 4, 2004,
266           pages 439-455.
267           <http://personal.strath.ac.uk/sergey.kitaev/publications.html>,
268           <http://www.jalc.de/issues/2004/issue_9_4/jalc-2004-439-455.php>
269
270           (Preprint as Sergey Kitaev and Toufik Mansour, "The Peano Curve and
271           Counting Occurrences of Some Patterns", October 2002.
272           <http://arxiv.org/abs/math/0210268/>.)
273
274       Their form is based on keeping the top-most U shape fixed and expanding
275       sub-parts.  This means the end segments alternate vertical and
276       horizontal in successive expansion levels.
277
278           direction            k=1              2       2
279             1 to 4                            *---*   *---*
280                                  2           1|  3|   |1  |3
281               1                *---*          *   *---*   *
282               |               1|   |3        1| 4   2   4 |3
283           4--- ---2            *   *          *---*   *---*
284               |                                  1|   |3       k=2
285               3                               *---*   *---*
286                                                 2       2
287
288           count segments in direction, for k >= 1
289           d(1,k) = 4^(k-1)                = 1,4,16,64,256,1024,4096,...
290           d(2,k) = 4^(k-1) + 2^(k-1) - 1  = 1,5,19,71,271,1055,4159,...
291           d(3,k) = 4^(k-1)                = 1,4,16,64,256,1024,4096,...
292           d(4,k) = 4^(k-1) - 2^(k-1)      = 0,2,12,56,240, 992,4032,...
293                                    (A000302, A099393, A000302, A020522)
294
295           total segments d(1,k)+d(2,k)+d(3,k)+d(4,k) = 4^k - 1
296
297       The path form here keeps the first segment direction fixed.  This means
298       a transpose 1<->2 and 3<->4 in odd levels.  The result is to take the
299       alternate d values as follows.  For k=0 there is a single point N=0 so
300       no line segments at all and so c(dir,0)=0.
301
302           first 4^k-1 segments
303
304           c(1,k) = / 0                        if k=0
305            North   | 4^(k-1) + 2^(k-1) - 1    if k odd >= 1
306                    \ 4^(k-1)                  if k even >= 2
307             = 0, 1, 4, 19, 64, 271, 1024, 4159, 16384, ...
308
309
310           c(2,k) = / 0                        if k=0
311            East    | 4^(k-1)                  if k odd >= 1
312                    \ 4^(k-1) + 2^(k-1) - 1    if k even >= 2
313             = 0, 1, 5, 16, 71, 256, 1055, 4096, 16511, ...
314
315           c(3,k) = / 0                        if k=0
316            South   | 4^(k-1) - 2^(k-1)        if k odd >= 1
317                    \ 4^(k-1)                  if k even >= 2
318             = 0, 0, 4, 12, 64, 240, 1024, 4032, 16384, ...
319
320           c(4,k) = / 0                        if k=0
321            West    | 4^(k-1)                  if k odd >= 1
322                    \ 4^(k-1) - 2^(k-1)        if k even >= 2
323             = 0, 1, 2, 16, 56, 256, 992, 4096, 16256, ...
324
325       The segment N=4^k-1 to N=4^k is North (direction 1) when k odd, or East
326       (direction 2) when k even.  That could be added to the respective cases
327       in c(1,k) and c(2,k) if desired.
328
329   Hamming Distance
330       The Hamming distance between integers X and Y is the number of bit
331       positions where the two values differ when written in binary.  On the
332       Hilbert curve each bit-pair of N becomes a bit of X and a bit of Y,
333
334              N      X   Y
335           ------   --- ---
336           0 = 00    0   0
337           1 = 01    1   0     <- difference 1 bit
338           2 = 10    1   1
339           3 = 11    0   1     <- difference 1 bit
340
341       So the Hamming distance for N=0to3 is 1 at N=1 and N=3.  At higher
342       levels, these X,Y bits may be transposed (swapped) or rotated by 180 or
343       both.  A transpose swapping X<->Y doesn't change the bit difference.  A
344       rotate by 180 is a flip 0<->1 of the bit in both X and Y, so that
345       doesn't change the bit difference either.
346
347       On that basis, the Hamming distance X,Y is the number of base4 digits
348       of N which are 01 or 11.  If bit positions are counted from 0 for the
349       least significant bit then
350
351           X,Y coordinates of N
352           HammingDist(X,Y) = count 1-bits at even bit positions in N
353                  = 0,1,0,1, 1,2,1,2, 0,1,0,1, 1,2,1,2, ... (A139351)
354
355       See also "Hamming Distance" in Math::PlanePath::CornerReplicate which
356       is the same formula, but arising directly from 01 or 11, no transpose
357       or rotate.
358

OEIS

360       This path is in Sloane's OEIS in several forms,
361
362           <http://oeis.org/A059252> (etc)
363
364           A059253    X coord
365           A059252    Y coord
366           A059261    X+Y
367           A059285    X-Y
368           A163547    X^2+Y^2 = radius squared
369           A139351    HammingDist(X,Y)
370           A059905    X xor Y, being ZOrderCurve X
371
372           A163365    sum N on diagonal
373           A163477    sum N on diagonal, divided by 4
374           A163482    N values on X axis
375           A163483    N values on Y axis
376           A062880    N values on diagonal X=Y (digits 0,2 in base 4)
377
378           A163538    dX -1,0,1 change in X
379           A163539    dY -1,0,1 change in Y
380           A163540    absolute direction of each step (0=E,1=N,2=W,3=S)
381           A163541    absolute direction, swapped X,Y
382           A163542    relative direction (ahead=0,right=1,left=2)
383           A163543    relative direction, swapped X,Y
384
385           A083885    count East segments N=0 to N=4^k (first 4^k segs)
386
387           A163900    distance dX^2+dY^2 between Hilbert and ZOrder
388           A165464    distance dX^2+dY^2 between Hilbert and Peano
389           A165466    distance dX^2+dY^2 between Hilbert and transposed Peano
390           A165465    N where Hilbert and Peano have same X,Y
391           A165467    N where Hilbert and Peano have transposed same X,Y
392
393       The following take points of the plane in various orders, each value in
394       the sequence being the N of the Hilbert curve at those positions.
395
396           A163355    N by the ZOrderCurve points sequence
397           A163356      inverse, ZOrderCurve by Hilbert points order
398           A166041    N by the PeanoCurve points sequence
399           A166042      inverse, PeanoCurve N by Hilbert points order
400           A163357    N by diagonals like Math::PlanePath::Diagonals with
401                        first Hilbert step along same axis the diagonals start
402           A163358      inverse
403           A163359    N by diagonals, transposed start along the opposite axis
404           A163360      inverse
405           A163361    A163357 + 1, numbering the Hilbert N's from N=1
406           A163362      inverse
407           A163363    A163355 + 1, numbering the Hilbert N's from N=1
408           A163364      inverse
409
410       The above sequences are permutations of the integers since all X,Y
411       positions of the first quadrant are covered by each path (Hilbert,
412       ZOrder, Peano).  The inverse permutations can be thought of taking X,Y
413       positions in the Hilbert order and asking what N the ZOrder, Peano or
414       Diagonals path would put there.
415
416       The A163355 permutation by ZOrderCurve can be considered for repeats or
417       cycles,
418
419           A163905    ZOrderCurve permutation A163355 applied twice
420           A163915    ZOrderCurve permutation A163355 applied three times
421           A163901    fixed points (N where X,Y same in both curves)
422           A163902    2-cycle points
423           A163903    3-cycle points
424           A163890    cycle lengths, points by N
425           A163904    cycle lengths, points by diagonals
426           A163910    count of cycles in 4^k blocks
427           A163911    max cycle length in 4^k blocks
428           A163912    LCM of cycle lengths in 4^k blocks
429           A163914    count of 3-cycles in 4^k blocks
430           A163909      those counts for even k only
431           A163891    N of previously unseen cycle length
432           A163893      first differences of those A163891
433           A163894    smallest value not an n-cycle
434           A163895      position of new high in A163894
435           A163896      value of new high in A163894
436
437           A163907    ZOrderCurve permutation twice, on points by diagonals
438           A163908      inverse of this
439
440       See examples/hilbert-oeis.pl for a sample program printing the A163359
441       permutation values.
442

SEE ALSO

444       Math::PlanePath, Math::PlanePath::HilbertSides,
445       Math::PlanePath::HilbertSpiral
446
447       Math::PlanePath::PeanoCurve, Math::PlanePath::ZOrderCurve,
448       Math::PlanePath::BetaOmega, Math::PlanePath::KochCurve
449
450       Math::Curve::Hilbert, Algorithm::SpatialIndex::Strategy::QuadTree
451

HOME PAGE

453       <http://user42.tuxfamily.org/math-planepath/index.html>
454

LICENSE

456       Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019,
457       2020, 2021 Kevin Ryde
458
459       This file is part of Math-PlanePath.
460
461       Math-PlanePath is free software; you can redistribute it and/or modify
462       it under the terms of the GNU General Public License as published by
463       the Free Software Foundation; either version 3, or (at your option) any
464       later version.
465
466       Math-PlanePath is distributed in the hope that it will be useful, but
467       WITHOUT ANY WARRANTY; without even the implied warranty of
468       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
469       General Public License for more details.
470
471       You should have received a copy of the GNU General Public License along
472       with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.
473
474
475
476perl v5.38.0                      2023-07-20  Math::PlanePath::HilbertCurve(3)
Impressum