File Coverage

blib/lib/Math/PlanePath/HilbertSpiral.pm
Criterion Covered Total %
statement 58 140 41.4
branch 6 58 10.3
condition n/a
subroutine 14 16 87.5
pod 3 3 100.0
total 81 217 37.3


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
2              
3             # This file is part of Math-PlanePath.
4             #
5             # Math-PlanePath is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-PlanePath is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-PlanePath. If not, see .
17              
18              
19             package Math::PlanePath::HilbertSpiral;
20 1     1   1264 use 5.004;
  1         3  
21 1     1   5 use strict;
  1         1  
  1         35  
22              
23 1     1   6 use vars '$VERSION', '@ISA';
  1         1  
  1         54  
24             $VERSION = 127;
25 1     1   607 use Math::PlanePath;
  1         3  
  1         26  
26 1     1   437 use Math::PlanePath::Base::NSEW;
  1         2  
  1         34  
27             @ISA = ('Math::PlanePath::Base::NSEW',
28             'Math::PlanePath');
29              
30             use Math::PlanePath::Base::Generic
31 1         38 'is_infinite',
32 1     1   6 'round_nearest';
  1         2  
33             use Math::PlanePath::Base::Digits
34 1     1   456 'digit_split_lowtohigh';
  1         2  
  1         48  
35              
36 1     1   480 use Math::PlanePath::BetaOmega 52;
  1         15  
  1         36  
37             *_y_round_down_len_level = \&Math::PlanePath::BetaOmega::_y_round_down_len_level;
38              
39             # uncomment this to run the ### lines
40             #use Smart::Comments;
41              
42              
43 1     1   5 use constant n_start => 0;
  1         2  
  1         45  
44 1     1   4 use constant xy_is_visited => 1;
  1         3  
  1         35  
45 1     1   4 use constant x_negative_at_n => 4;
  1         2  
  1         43  
46 1     1   5 use constant y_negative_at_n => 8;
  1         2  
  1         862  
47              
48              
49             #------------------------------------------------------------------------------
50              
51             # generated by tools/hilbert-spiral-table.pl
52             #
53             my @next_state = (8,0,0,12, 12,4,4,8, 0,8,8,4, 4,12,12,0,
54             20,0,0,12, 16,4,4,8);
55             my @digit_to_x = (0,1,1,0, 1,0,0,1, 0,0,1,1, 1,1,0,0,
56             0,1,1,0, 1,0,0,1);
57             my @digit_to_y = (0,0,1,1, 1,1,0,0, 0,1,1,0, 1,0,0,1,
58             0,0,1,1, 1,1,0,0);
59             my @xy_to_digit = (0,3,1,2, 2,1,3,0, 0,1,3,2, 2,3,1,0,
60             0,3,1,2, 2,1,3,0);
61             my @min_digit = (0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
62             2,2,3,1, 0,0,1,0, 0,undef,undef,undef,
63             0,0,3,0, 0,2,1,1, 2,undef,undef,undef,
64             2,1,1,2, 0,0,3,0, 0,undef,undef,undef,
65             0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
66             2,2,3,1, 0,0,1,0, 0);
67             my @max_digit = (0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
68             2,3,3,2, 3,3,1,1, 0,undef,undef,undef,
69             0,3,3,1, 3,3,1,2, 2,undef,undef,undef,
70             2,2,1,3, 3,1,3,3, 0,undef,undef,undef,
71             0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
72             2,3,3,2, 3,3,1,1, 0);
73             # neg state 20
74              
75             sub n_to_xy {
76 2161     2161 1 18932 my ($self, $n) = @_;
77             ### HilbertSpiral n_to_xy(): $n
78             ### hex: sprintf "%#X", $n
79              
80 2161 50       3512 if ($n < 0) { return; }
  0         0  
81 2161 50       3465 if (is_infinite($n)) { return ($n,$n); }
  0         0  
82              
83 2161         3262 my $int = int($n);
84 2161         2465 $n -= $int;
85              
86 2161         3272 my @digits = digit_split_lowtohigh($int,4);
87 2161         3696 my $len = ($n*0 + 2) ** scalar(@digits); # inherit possible bigint 1
88              
89 2161 100       3254 my $state = ($#digits & 1 ? 4 : 0);
90 2161         2408 my $dir = $state + 2; # default if all $digit==3
91             ### @digits
92              
93 2161         2475 my $x = my $y = 0;
94              
95 2161         3597 while (defined (my $digit = pop @digits)) { # high to low
96 10091         11302 $len /= 2;
97 10091         10695 $state += $digit;
98 10091 100       13290 if ($digit != 3) {
99 7410         7964 $dir = $state; # lowest non-3 digit
100             }
101              
102             ### at: "$x,$y len=$len"
103             ### $state
104             ### $dir
105             ### digit_to_x: $digit_to_x[$state]
106             ### digit_to_y: $digit_to_y[$state]
107             ### next_state: $next_state[$state]
108              
109 10091         11739 my $offset = scalar(@digits) & 1;
110 10091         12747 $x += $len * ($digit_to_x[$state] - $offset);
111 10091         11254 $y += $len * ($digit_to_y[$state] - $offset);
112 10091         16026 $state = $next_state[$state];
113             }
114              
115              
116             ### frac: $n
117             ### $dir
118             ### dir dx: ($digit_to_x[$dir+1] - $digit_to_x[$dir])
119             ### dir dy: ($digit_to_y[$dir+1] - $digit_to_y[$dir])
120             ### x: $n * ($digit_to_x[$dir+1] - $digit_to_x[$dir]) + $x
121             ### y: $n * ($digit_to_y[$dir+1] - $digit_to_y[$dir]) + $y
122              
123             # with $n fractional part
124 2161         5742 return ($n * ($digit_to_x[$dir+1] - $digit_to_x[$dir]) + $x,
125             $n * ($digit_to_y[$dir+1] - $digit_to_y[$dir]) + $y);
126             }
127              
128             sub xy_to_n {
129 0     0 1   my ($self, $x, $y) = @_;
130             ### HilbertSpiral xy_to_n(): "$x, $y"
131              
132 0           $x = round_nearest ($x);
133 0           $y = round_nearest ($y);
134              
135 0           my $n = ($x * 0 * $y);
136              
137 0           my ($len, $level) = _y_round_down_len_level ($x);
138             {
139 0           my ($ylen, $ylevel) = _y_round_down_len_level ($y);
  0            
140             ### y len/level: "$ylen $ylevel"
141 0 0         if ($ylevel > $level) {
142 0           $level = $ylevel;
143 0           $len = $ylen;
144             }
145             }
146 0 0         if (is_infinite($len)) {
147 0           return $len;
148             }
149              
150             ### $len
151             ### $level
152              
153 0           my $state;
154             {
155 0           my $offset;
  0            
156 0 0         if ($level & 1) {
157 0           $state = 4;
158 0           $offset = 4*$len;
159             } else {
160 0           $state = 0;
161 0           $offset = 2*$len;
162             }
163 0           $offset -= 2;
164 0           $offset /= 3;
165 0           $y += $offset;
166 0           $x += $offset;
167             # $x,$y now relative to Xmin(level),Ymin(level),
168             # so in range 0 <= $x,$y < 2*len
169             }
170             ### offset x,y to: "$x, $y"
171              
172 0           for (;;) {
173             ### at: "$x,$y len=$len"
174             ### assert: $x >= 0
175             ### assert: $y >= 0
176             ### assert: $x < 2*$len
177             ### assert: $y < 2*$len
178              
179 0           my $xo;
180 0 0         if ($xo = ($x >= $len)) {
181 0           $x -= $len;
182             }
183 0           my $yo;
184 0 0         if ($yo = ($y >= $len)) {
185 0           $y -= $len;
186             }
187             ### xy bits: ($xo+0).", ".($yo+0)
188              
189 0           my $digit = $xy_to_digit[$state + 2*$xo + $yo];
190 0           $n = 4*$n + $digit;
191 0           $state = $next_state[$state+$digit];
192              
193 0 0         last if --$level < 0;
194 0           $len /= 2;
195             }
196              
197             ### assert: $x == 0
198             ### assert: $y == 0
199              
200 0           return $n;
201             }
202              
203              
204             # This finds the exact minimum/maximum N in the given rectangle.
205             #
206             # The strategy is similar to xy_to_n(), except that at each bit position
207             # instead of taking a bit of x,y from the input instead those bits are
208             # chosen from among the 4 sub-parts according to which has the maximum N and
209             # is within the given target rectangle. The final result is both an $n_max
210             # and a $x_max,$y_max which is its position, but only the $n_max is
211             # returned.
212             #
213             # At a given sub-part the comparisons ask whether x1 is above or below the
214             # midpoint, and likewise x2,y1,y2. Since x2>=x1 and y2>=y1 there's only 3
215             # combinations of x1>=cmp,x2>=cmp, not 4.
216              
217             # exact
218             sub rect_to_n_range {
219 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
220             ### HilbertSpiral rect_to_n_range(): "$x1,$y1, $x2,$y2"
221              
222 0           $x1 = round_nearest ($x1);
223 0           $y1 = round_nearest ($y1);
224 0           $x2 = round_nearest ($x2);
225 0           $y2 = round_nearest ($y2);
226 0 0         ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
227 0 0         ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
228              
229             # If y1/y2 both positive or both negative then only look at the bigger of
230             # the two. If y1 negative and y2 positive then consider both.
231 0           my $len = 1;
232 0           my $level = 0;
233 0 0         foreach my $z (($x2 > 0 ? ($x2) : ()),
    0          
    0          
    0          
234             ($x1 < 0 ? ($x1) : ()),
235             ($y2 > 0 ? ($y2) : ()),
236             ($y1 < 0 ? ($y1) : ())) {
237 0           my ($zlen, $zlevel) = _y_round_down_len_level ($z);
238             ### y len/level: "$zlen $zlevel"
239 0 0         if ($zlevel > $level) {
240 0           $level = $zlevel;
241 0           $len = $zlen;
242             }
243             }
244 0 0         if (is_infinite($len)) {
245 0           return (0, $len);
246             }
247              
248             # At this point an easy over-estimate would be:
249             # return (0, $len*$len*4-1);
250              
251 0           my $n_min = my $n_max = 0;
252 0           my $x_min = my $x_max = my $y_min = my $y_max
253             = - (4**int(($level+1)/2) - 1) * 2 / 3;
254 0 0         my $min_state = my $max_state = ($level & 1 ? 20 : 16);
255             ### $x_min
256             ### $y_min
257              
258 0           while ($level >= 0) {
259             ### $level
260             ### $len
261             {
262 0           my $x_cmp = $x_min + $len;
263 0           my $y_cmp = $y_min + $len;
264 0 0         my $digit = $min_digit[3*$min_state
    0          
    0          
    0          
265             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
266             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];
267              
268 0           $n_min = 4*$n_min + $digit;
269 0           $min_state += $digit;
270 0 0         if ($digit_to_x[$min_state]) { $x_min += $len; }
  0            
271 0           $y_min += $len * $digit_to_y[$min_state];
272 0           $min_state = $next_state[$min_state];
273             }
274             {
275 0           my $x_cmp = $x_max + $len;
  0            
  0            
276 0           my $y_cmp = $y_max + $len;
277 0 0         my $digit = $max_digit[3*$max_state
    0          
    0          
    0          
278             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
279             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];
280              
281 0           $n_max = 4*$n_max + $digit;
282 0           $max_state += $digit;
283 0 0         if ($digit_to_x[$max_state]) { $x_max += $len; }
  0            
284 0           $y_max += $len * $digit_to_y[$max_state];
285 0           $max_state = $next_state[$max_state];
286             }
287              
288 0           $len = int($len/2);
289 0           $level--;
290             }
291              
292 0           return ($n_min, $n_max);
293             }
294              
295             #------------------------------------------------------------------------------
296             # levels
297              
298 1     1   6 use Math::PlanePath::HilbertCurve;
  1         2  
  1         55  
299             *level_to_n_range = \&Math::PlanePath::HilbertCurve::level_to_n_range;
300             *n_to_level = \&Math::PlanePath::HilbertCurve::n_to_level;
301              
302             #------------------------------------------------------------------------------
303             1;
304             __END__