1Math::PlanePath::KochCuUrsveer(3C)ontributed Perl DocumeMnattaht:i:oPnlanePath::KochCurve(3)
2
6 Math::PlanePath::KochCurve -- horizontal Koch curve
7
9 use Math::PlanePath::KochCurve;
10 my $path = Math::PlanePath::KochCurve->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This is an integer version of the self-similar Koch curve,
15 Helge von Koch, "Une Méthode Géométrique Élémentaire pour l'Étude
17 de Certaines Questions de la Théorie des Courbes Planes", Acta
18 Arithmetica, volume 30, 1906, pages 145-174.
19 <http://archive.org/details/actamathematica11lefgoog>
20 It goes along the X axis and makes triangular excursions upwards.
22 8 3
24 / \
25 6---- 7 9----10 18-... 2
26 \ / \
27 2 5 11 14 17 1
28 / \ / \ / \ /
29 0----1 3---- 4 12----13 15----16 <- Y=0
30 ^
32 X=0 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
33 The replicating shape is the initial N=0 to N=4,
35 *
37 / \
38 *---* *---*
39 which is rotated and repeated 3 times in the same pattern to give
41 sections N=4 to N=8, N=8 to N=12, and N=12 to N=16. Then that N=0 to
42 N=16 is itself replicated three times at the angles of the base
43 pattern, and so on infinitely.
44 The X,Y coordinates are arranged on a square grid using every second
46 point, per "Triangular Lattice" in Math::PlanePath. The result is
47 flattened triangular segments with diagonals at a 45 degree angle.
48 Level Ranges
50 Each replication adds 3 copies of the existing points and is thus 4
51 times bigger, so if N=0 to N=4 is reckoned as level 1 then a given
52 replication level goes from
53 Nstart = 0
55 Nlevel = 4^level (inclusive)
56 Each replication is 3 times the width. The initial N=0 to N=4 figure
58 is 6 wide and in general a level runs from
59 Xstart = 0
61 Xlevel = 2*3^level at N=Nlevel
62 The highest Y is 3 times greater at each level similarly. The peak is
64 at the midpoint of each level,
65 Npeak = (4^level)/2
67 Ypeak = 3^level
68 Xpeak = 3^level
69 It can be seen that the N=6 point backtracks horizontally to the same X
71 as the start of its section N=4 to N=8. This happens in the further
72 replications too and is the maximum extent of the backtracking.
73 The Nlevel is multiplied by 4 to get the end of the next higher level.
75 The same 4*N can be applied to all points N=0 to N=Nlevel to get the
76 same shape but a factor of 3 bigger X,Y coordinates. The in-between
77 points 4*N+1, 4*N+2 and 4*N+3 are then new finer structure in the
78 higher level.
79 Fractal
81 Koch conceived the curve as having a fixed length and infinitely fine
82 structure, making it continuous everywhere but differentiable nowhere.
83 The code here can be pressed into use for that sort of construction for
84 a given level of granularity by scaling
85 X/3^level
87 Y/3^level
88 which makes it a fixed 2 wide by 1 high. Or for unit-side equilateral
90 triangles then apply further factors 1/2 and sqrt(3)/2, as noted in
91 "Triangular Lattice" in Math::PlanePath.
92 (X/2) / 3^level
94 (Y*sqrt(3)/2) / 3^level
95 Area
97 The area under the curve to a given level can be calculated from its
98 self-similar nature. The curve at level+1 is 3 times wider and higher
99 and adds a triangle of unit area onto each line segment. So reckoning
100 the line segment N=0 to N=1 as level=0 (which is area[0]=0),
101 area[level] = 9*area[level-1] + 4^(level-1)
103 = 4^(level-1) + 9*4^(level-2) + ... + 9^(level-1)*4^0
104 9^level - 4^level
106 = -----------------
107 5
108 = 0, 1, 13, 133, 1261, 11605, 105469, ... (A016153)
110 The sides are 6 different angles. The triangles added on the sides are
112 always the same shape either pointing up or down. Base width=2 and
113 height=1 gives area=1.
114 * *-----* ^
116 / \ \ / | height=1
117 / \ \ / |
118 *-----* * v triangle area = 2*1/2 = 1
119 <-----> width=2
121 If the Y coordinates are stretched to make equilateral triangles then
123 the number of triangles is not changed and so the area increases by a
124 factor of the area of the equilateral triangle, sqrt(3)/4.
125
127 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
128 classes.
129 "$path = Math::PlanePath::KochCurve->new ()"
131 Create and return a new path object.
132 "($x,$y) = $path->n_to_xy ($n)"
134 Return the X,Y coordinates of point number $n on the path. Points
135 begin at 0 and if "$n < 0" then the return is an empty list.
136 Fractional positions give an X,Y position along a straight line
138 between the integer positions.
139 "($n_lo, $n_hi) = $path->rect_to_n_range ($x1,$y1, $x2,$y2)"
141 The returned range is exact, meaning $n_lo and $n_hi are the
142 smallest and biggest in the rectangle.
143 "$n = $path->n_start()"
145 Return 0, the first N in the path.
146 Level Methods
148 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
149 Return "(0, 4**$level)".
150
152 N to Turn
153 The curve always turns either +60 degrees or -120 degrees, it never
154 goes straight ahead. In the base 4 representation of N, the lowest
155 non-zero digit gives the turn. The first turn is at N=1 so there's
156 always a non-zero digit in N.
157 low digit
159 base 4 turn
160 --------- ------------
161 1 +60 degrees (left)
162 2 -120 degrees (right)
163 3 +60 degrees (left)
164 For example N=8 is 20 base 4, so lowest nonzero "2" means turn -120
166 degrees for the next segment.
167 If the least significant digit is non-zero then it determines the turn,
169 making the base N=0 to N=4 shape. If the least significant is zero
170 then the next level up is in control, eg. N=0,4,8,12,16, making a turn
171 according to the base shape again at that higher level. The first and
172 last segments of the base shape are "straight" so there's no extra
173 adjustment to apply in those higher digits.
174 This base 4 digit rule is equivalent to counting low 0-bits. A low
176 base-4 digit 1 or 3 is an even number of low 0-bits and a low digit 2
177 is an odd number of low 0-bits.
178 count low 0-bits turn
180 ---------------- ------------
181 even +60 degrees (left)
182 odd -120 degrees (right)
183 For example N=8 in binary "1000" has 3 low 0-bits and 3 is odd so turn
185 -120 degrees (right).
186 See "Turn" in Math::PlanePath::GrayCode for a similar turn sequence
188 arising from binary Gray code.
189 N to Next Turn
191 The turn at N+1, ie the next turn, can be found from the base-4 digits
192 by considering how the digits of N change when 1 is added, and the low-
193 digit turn calculation is applied on those changed digits.
194 Adding 1 means low digit 0, 1 or 2 will become non-zero. Any low 3s
196 wrap around to become low 0s. So the turn at N+1 can be found from the
197 digits of N by seeking the lowest non-3
198 lowest non-3 turn
200 digit of N at N+1
201 ------------ ------------
202 0 +60 degrees (left)
203 1 -120 degrees (right)
204 2 +60 degrees (left)
205 N to Direction
207 The total turn at a given N can be found by counting digits 1 and 2 in
208 base 4.
209 direction = ((count of 1-digits in base 4)
211 - (count of 2-digits in base 4)) * 60 degrees
212 For example N=11 is "23" in base 4, so 60*(0-1) = -60 degrees.
214 In this formula the count of 1s and 2s can go past 360 degrees,
216 representing a spiralling around which occurs at progressively higher
217 replication levels. The direction can be taken mod 360 degrees, or the
218 count mod 6, for a direction 0 to 5 if desired.
219 N to abs(dX),abs(dY)
221 The direction expressed as abs(dX) and abs(dY) can be calculated simply
222 from N modulo 3. abs(dX) is a repeating pattern 2,1,1 and abs(dY)
223 repeating 0,1,1.
224 N mod 3 abs(dX),abs(dY)
226 ------- ---------------
227 0 2,0 horizontal, East or West
228 1 1,1 slope North-East or South-West
229 2 1,1 slope North-West or South-East
230 This works because the direction calculation above corresponds to N mod
232 3. Each N digit in base 4 becomes
233 N digit
235 base 4 direction add
236 ------- -------------
237 0 0
238 1 1
239 2 -1
240 3 0
241 Notice that direction == Ndigit mod 3. Then because 4==1 mod 3 the
243 power-of-4 for each digit reduces down to 1,
244 N = 4^k * digit_k + ... 4^0 * digit_0
246 N mod 3 = 1 * digit_k + ... 1 * digit_0
247 = digit_k + ... digit_0
248 same as
249 direction = digit_k + ... + digit_0 taken mod 3
250 Rectangle to N Range -- Level
252 An easy over-estimate of the N values in a rectangle can be had from
253 the Xlevel formula above. If Xlevel>rectangleX then Nlevel is past the
254 rectangle extent.
255 X = 2*3^level
257 so
259 floorlevel = floor log_base_3(X/2)
261 Nhi = 4^(floorlevel+1) - 1
262 For example a rectangle extending to X=13 has floorlevel =
264 floor(log3(13/2))=1 and so Nhi=4^(1+1)-1=15.
265 The rounding-down of the log3 ensures a point such as X=18 which is the
267 first in the next Nlevel will give that next level. So
268 floorlevel=log3(18/2)=2 (exactly) and Nhi=4^(2+1)-1=63.
269 The worst case for this over-estimate is when rectangleX==Xlevel, ie.
271 the rectangle is just into the next level. In that case Nhi is nearly
272 a factor 4 bigger than it needs to be.
273 Rectangle to N Range -- Exact
275 The exact Nlo and Nhi in a rectangle can be found by searching along
276 the curve. For Nlo search forward from the origin N=0. For Nhi search
277 backward from the Nlevel over-estimate described above.
278 At a given digit position in the prospective N the sub-part of the
280 curve comprising the lower digits has a certain triangular extent. If
281 it's outside the target rectangle then step to the next digit value,
282 and to the next of the digit above when past digit=3 (or below digit=0
283 when searching backwards).
284 There's six possible orientations for the curve sub-part. In the
286 following pictures "o" is the start and the surrounding lines show the
287 triangular extent. There's just four curve parts shown in each, but
288 these triangles bound a sub-curve of any level.
289 rot=0 -+- +-----------------+
291 -- -- - .-+-* *-+-o -
292 -- * -- -- \ / --
293 -- / \ -- -- * --
294 - o-+-* *-+-. - -- --
295 +-----------------+ rot=3 -+-
296 rot=1
298 +---------+ rot=4 /+
299 | . / / |
300 | / / / o|
301 |*-+-* / / / |
302 | \ / / * |
303 | * / / \ |
304 | / / / *-+-*|
305 |o / / / |
306 | / / . |
307 +/ +---------+
308 +\ rot=2 +---------+
310 | \ \ o |
311 |. \ \ \ |
312 | \ \ \ *-+-*|
313 | * \ \ / |
314 | / \ \ * |
315 |*-+-* \ \ \ |
316 | \ \ \ .|
317 | o \ rot=5 \ |
318 +---------+ \+
319 The "." is the start of the next sub-curve. It belongs to the next
321 digit value and so can be excluded. For rot=0 and rot=3 this means
322 simply shortening the X range permitted. For rot=1 and rot=4 similarly
323 the Y range. For rot=2 and rot=5 it would require a separate test.
324 Tight sub-part extent checking reduces the sub-parts which are
326 examined, but it works perfectly well with a looser check, such as a
327 square box for the sub-curve extents. Doing that might be easier if
328 the target region is not a rectangle but instead some trickier shape.
329
331 The Koch curve is in Sloane's Online Encyclopedia of Integer Sequences
332 in various forms,
333 <http://oeis.org/A035263> (etc)
335 A335358 (X-Y)/2 diagonal coordinate
337 A335359 Y coordinate
338 A035263 turn 1=left,0=right, by morphism
340 A096268 turn 0=left,1=right, period doubling sequence
341 A056832 turn 1=left,2=right, by replicate and flip last
342 A309873 turn 1=left,-1=right
343 A029883 turn +/-1=left,0=right, Thue-Morse first differences
344 A089045 turn +/-1=left,0=right, by +/- something
345 A177702 abs(dX) from N=1 onwards, being 1,1,2 repeating
347 A011655 abs(dY), being 0,1,1 repeating
348 A003159 N positions of left turns, ending even number 0 bits
350 A036554 N positions of right turns, ending odd number 0 bits
351 A332206 N on X axis
352 A001196 N segments on X axis (N and N+1 on X axis)
353 A065359 segment direction, *60 degrees
355 A229216 segment direction, 1,2,3,-1,-2,-3
356 A050292 num left turns 1 to N
357 A123087 num right turns 1 to N
358 A020988 num left turns 1 to 4^k-1, being 2*(4^k-1)/3
359 A002450 num right turns 1 to 4^k-1, being (4^k-1)/3
360 A016153 area under the curve, (9^k-4^k)/5
361 For reference, A217586 is not quite the same as A096268 right turn.
363 A217586 differs by a 0<->1 flip at N=2^k due to different initial
364 a(1)=1.
365
367 Math::PlanePath, Math::PlanePath::PeanoCurve,
368 Math::PlanePath::HilbertCurve, Math::PlanePath::KochPeaks,
369 Math::PlanePath::KochSnowflakes, Math::PlanePath::CCurve
370 Math::Fractal::Curve
372
374 <http://user42.tuxfamily.org/math-planepath/index.html>
375
377 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
378 Kevin Ryde
379 Math-PlanePath is free software; you can redistribute it and/or modify
381 it under the terms of the GNU General Public License as published by
382 the Free Software Foundation; either version 3, or (at your option) any
383 later version.
384 Math-PlanePath is distributed in the hope that it will be useful, but
386 WITHOUT ANY WARRANTY; without even the implied warranty of
387 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
388 General Public License for more details.
389 You should have received a copy of the GNU General Public License along
391 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
392perl v5.36.0 2022-07-22 Math::PlanePath::KochCurve(3)