cebf0b0d |
1 | /* penrose.c |
2 | * |
3 | * Penrose tile generator. |
4 | * |
5 | * Uses half-tile technique outlined on: |
6 | * |
7 | * http://tartarus.org/simon/20110412-penrose/penrose.xhtml |
8 | */ |
9 | |
10 | #include <assert.h> |
11 | #include <string.h> |
12 | #include <math.h> |
13 | #include <stdio.h> |
14 | |
15 | #include "puzzles.h" /* for malloc routines, and PI */ |
16 | #include "penrose.h" |
17 | |
18 | /* ------------------------------------------------------- |
19 | * 36-degree basis vector arithmetic routines. |
20 | */ |
21 | |
22 | /* Imagine drawing a |
23 | * ten-point 'clock face' like this: |
24 | * |
25 | * -E |
26 | * -D | A |
27 | * \ | / |
28 | * -C. \ | / ,B |
29 | * `-._\|/_,-' |
30 | * ,-' /|\ `-. |
31 | * -B' / | \ `C |
32 | * / | \ |
33 | * -A | D |
34 | * E |
35 | * |
36 | * In case the ASCII art isn't clear, those are supposed to be ten |
37 | * vectors of length 1, all sticking out from the origin at equal |
38 | * angular spacing (hence 36 degrees). Our basis vectors are A,B,C,D (I |
39 | * choose them to be symmetric about the x-axis so that the final |
40 | * translation into 2d coordinates will also be symmetric, which I |
41 | * think will avoid minor rounding uglinesses), so our vector |
42 | * representation sets |
43 | * |
44 | * A = (1,0,0,0) |
45 | * B = (0,1,0,0) |
46 | * C = (0,0,1,0) |
47 | * D = (0,0,0,1) |
48 | * |
49 | * The fifth vector E looks at first glance as if it needs to be |
50 | * another basis vector, but in fact it doesn't, because it can be |
51 | * represented in terms of the other four. Imagine starting from the |
52 | * origin and following the path -A, +B, -C, +D: you'll find you've |
53 | * traced four sides of a pentagram, and ended up one E-vector away |
54 | * from the origin. So we have |
55 | * |
56 | * E = (-1,1,-1,1) |
57 | * |
58 | * This tells us that we can rotate any vector in this system by 36 |
59 | * degrees: if we start with a*A + b*B + c*C + d*D, we want to end up |
60 | * with a*B + b*C + c*D + d*E, and we substitute our identity for E to |
61 | * turn that into a*B + b*C + c*D + d*(-A+B-C+D). In other words, |
62 | * |
63 | * rotate_one_notch_clockwise(a,b,c,d) = (-d, d+a, -d+b, d+c) |
64 | * |
65 | * and you can verify for yourself that applying that operation |
66 | * repeatedly starting with (1,0,0,0) cycles round ten vectors and |
67 | * comes back to where it started. |
68 | * |
69 | * The other operation that may be required is to construct vectors |
70 | * with lengths that are multiples of phi. That can be done by |
71 | * observing that the vector C-B is parallel to E and has length 1/phi, |
72 | * and the vector D-A is parallel to E and has length phi. So this |
73 | * tells us that given any vector, we can construct one which points in |
74 | * the same direction and is 1/phi or phi times its length, like this: |
75 | * |
76 | * divide_by_phi(vector) = rotate(vector, 2) - rotate(vector, 3) |
77 | * multiply_by_phi(vector) = rotate(vector, 1) - rotate(vector, 4) |
78 | * |
79 | * where rotate(vector, n) means applying the above |
80 | * rotate_one_notch_clockwise primitive n times. Expanding out the |
81 | * applications of rotate gives the following direct representation in |
82 | * terms of the vector coordinates: |
83 | * |
84 | * divide_by_phi(a,b,c,d) = (b-d, c+d-b, a+b-c, c-a) |
85 | * multiply_by_phi(a,b,c,d) = (a+b-d, c+d, a+b, c+d-a) |
86 | * |
87 | * and you can verify for yourself that those two operations are |
88 | * inverses of each other (as you'd hope!). |
89 | * |
90 | * Having done all of this, testing for equality between two vectors is |
91 | * a trivial matter of comparing the four integer coordinates. (Which |
92 | * it _wouldn't_ have been if we'd kept E as a fifth basis vector, |
93 | * because then (-1,1,-1,1,0) and (0,0,0,0,1) would have had to be |
94 | * considered identical. So leaving E out is vital.) |
95 | */ |
96 | |
97 | struct vector { int a, b, c, d; }; |
98 | |
fd66a01d |
99 | static vector v_origin(void) |
cebf0b0d |
100 | { |
101 | vector v; |
102 | v.a = v.b = v.c = v.d = 0; |
103 | return v; |
104 | } |
105 | |
106 | /* We start with a unit vector of B: this means we can easily |
107 | * draw an isoceles triangle centred on the X axis. */ |
108 | #ifdef TEST_VECTORS |
109 | |
fd66a01d |
110 | static vector v_unit(void) |
cebf0b0d |
111 | { |
112 | vector v; |
113 | |
114 | v.b = 1; |
115 | v.a = v.c = v.d = 0; |
116 | return v; |
117 | } |
118 | |
119 | #endif |
120 | |
121 | #define COS54 0.5877852 |
122 | #define SIN54 0.8090169 |
123 | #define COS18 0.9510565 |
124 | #define SIN18 0.3090169 |
125 | |
126 | /* These two are a bit rough-and-ready for now. Note that B/C are |
127 | * 18 degrees from the x-axis, and A/D are 54 degrees. */ |
128 | double v_x(vector *vs, int i) |
129 | { |
130 | return (vs[i].a + vs[i].d) * COS54 + |
131 | (vs[i].b + vs[i].c) * COS18; |
132 | } |
133 | |
134 | double v_y(vector *vs, int i) |
135 | { |
136 | return (vs[i].a - vs[i].d) * SIN54 + |
137 | (vs[i].b - vs[i].c) * SIN18; |
138 | |
139 | } |
140 | |
141 | static vector v_trans(vector v, vector trans) |
142 | { |
143 | v.a += trans.a; |
144 | v.b += trans.b; |
145 | v.c += trans.c; |
146 | v.d += trans.d; |
147 | return v; |
148 | } |
149 | |
150 | static vector v_rotate_36(vector v) |
151 | { |
152 | vector vv; |
153 | vv.a = -v.d; |
154 | vv.b = v.d + v.a; |
155 | vv.c = -v.d + v.b; |
156 | vv.d = v.d + v.c; |
157 | return vv; |
158 | } |
159 | |
160 | static vector v_rotate(vector v, int ang) |
161 | { |
162 | int i; |
163 | |
164 | assert((ang % 36) == 0); |
165 | while (ang < 0) ang += 360; |
166 | ang = 360-ang; |
167 | for (i = 0; i < (ang/36); i++) |
168 | v = v_rotate_36(v); |
169 | return v; |
170 | } |
171 | |
172 | #ifdef TEST_VECTORS |
173 | |
174 | static vector v_scale(vector v, int sc) |
175 | { |
176 | v.a *= sc; |
177 | v.b *= sc; |
178 | v.c *= sc; |
179 | v.d *= sc; |
180 | return v; |
181 | } |
182 | |
183 | #endif |
184 | |
185 | static vector v_growphi(vector v) |
186 | { |
187 | vector vv; |
188 | vv.a = v.a + v.b - v.d; |
189 | vv.b = v.c + v.d; |
190 | vv.c = v.a + v.b; |
191 | vv.d = v.c + v.d - v.a; |
192 | return vv; |
193 | } |
194 | |
195 | static vector v_shrinkphi(vector v) |
196 | { |
197 | vector vv; |
198 | vv.a = v.b - v.d; |
199 | vv.b = v.c + v.d - v.b; |
200 | vv.c = v.a + v.b - v.c; |
201 | vv.d = v.c - v.a; |
202 | return vv; |
203 | } |
204 | |
205 | #ifdef TEST_VECTORS |
206 | |
207 | static const char *v_debug(vector v) |
208 | { |
209 | static char buf[255]; |
210 | sprintf(buf, |
211 | "(%d,%d,%d,%d)[%2.2f,%2.2f]", |
212 | v.a, v.b, v.c, v.d, v_x(&v,0), v_y(&v,0)); |
213 | return buf; |
214 | } |
215 | |
216 | #endif |
217 | |
218 | /* ------------------------------------------------------- |
219 | * Tiling routines. |
220 | */ |
221 | |
fd66a01d |
222 | static vector xform_coord(vector vo, int shrink, vector vtrans, int ang) |
cebf0b0d |
223 | { |
224 | if (shrink < 0) |
225 | vo = v_shrinkphi(vo); |
226 | else if (shrink > 0) |
227 | vo = v_growphi(vo); |
228 | |
229 | vo = v_rotate(vo, ang); |
230 | vo = v_trans(vo, vtrans); |
231 | |
232 | return vo; |
233 | } |
234 | |
235 | |
236 | #define XFORM(n,o,s,a) vs[(n)] = xform_coord(v_edge, (s), vs[(o)], (a)) |
237 | |
238 | static int penrose_p2_small(penrose_state *state, int depth, int flip, |
239 | vector v_orig, vector v_edge); |
240 | |
241 | static int penrose_p2_large(penrose_state *state, int depth, int flip, |
242 | vector v_orig, vector v_edge) |
243 | { |
244 | vector vv_orig, vv_edge; |
245 | |
246 | #ifdef DEBUG_PENROSE |
247 | { |
248 | vector vs[3]; |
249 | vs[0] = v_orig; |
250 | XFORM(1, 0, 0, 0); |
251 | XFORM(2, 0, 0, -36*flip); |
252 | |
253 | state->new_tile(state, vs, 3, depth); |
254 | } |
255 | #endif |
256 | |
257 | if (flip > 0) { |
258 | vector vs[4]; |
259 | |
260 | vs[0] = v_orig; |
261 | XFORM(1, 0, 0, -36); |
262 | XFORM(2, 0, 0, 0); |
263 | XFORM(3, 0, 0, 36); |
264 | |
265 | state->new_tile(state, vs, 4, depth); |
266 | } |
267 | if (depth >= state->max_depth) return 0; |
268 | |
269 | vv_orig = v_trans(v_orig, v_rotate(v_edge, -36*flip)); |
270 | vv_edge = v_rotate(v_edge, 108*flip); |
271 | |
272 | penrose_p2_small(state, depth+1, flip, |
273 | v_orig, v_shrinkphi(v_edge)); |
274 | |
275 | penrose_p2_large(state, depth+1, flip, |
276 | vv_orig, v_shrinkphi(vv_edge)); |
277 | penrose_p2_large(state, depth+1, -flip, |
278 | vv_orig, v_shrinkphi(vv_edge)); |
279 | |
280 | return 0; |
281 | } |
282 | |
283 | static int penrose_p2_small(penrose_state *state, int depth, int flip, |
284 | vector v_orig, vector v_edge) |
285 | { |
286 | vector vv_orig; |
287 | |
288 | #ifdef DEBUG_PENROSE |
289 | { |
290 | vector vs[3]; |
291 | vs[0] = v_orig; |
292 | XFORM(1, 0, 0, 0); |
293 | XFORM(2, 0, -1, -36*flip); |
294 | |
295 | state->new_tile(state, vs, 3, depth); |
296 | } |
297 | #endif |
298 | |
299 | if (flip > 0) { |
300 | vector vs[4]; |
301 | |
302 | vs[0] = v_orig; |
303 | XFORM(1, 0, 0, -72); |
304 | XFORM(2, 0, -1, -36); |
305 | XFORM(3, 0, 0, 0); |
306 | |
307 | state->new_tile(state, vs, 4, depth); |
308 | } |
309 | |
310 | if (depth >= state->max_depth) return 0; |
311 | |
312 | vv_orig = v_trans(v_orig, v_edge); |
313 | |
314 | penrose_p2_large(state, depth+1, -flip, |
315 | v_orig, v_shrinkphi(v_rotate(v_edge, -36*flip))); |
316 | |
317 | penrose_p2_small(state, depth+1, flip, |
318 | vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip))); |
319 | |
320 | return 0; |
321 | } |
322 | |
323 | static int penrose_p3_small(penrose_state *state, int depth, int flip, |
324 | vector v_orig, vector v_edge); |
325 | |
326 | static int penrose_p3_large(penrose_state *state, int depth, int flip, |
327 | vector v_orig, vector v_edge) |
328 | { |
329 | vector vv_orig; |
330 | |
331 | #ifdef DEBUG_PENROSE |
332 | { |
333 | vector vs[3]; |
334 | vs[0] = v_orig; |
335 | XFORM(1, 0, 1, 0); |
336 | XFORM(2, 0, 0, -36*flip); |
337 | |
338 | state->new_tile(state, vs, 3, depth); |
339 | } |
340 | #endif |
341 | |
342 | if (flip > 0) { |
343 | vector vs[4]; |
344 | |
345 | vs[0] = v_orig; |
346 | XFORM(1, 0, 0, -36); |
347 | XFORM(2, 0, 1, 0); |
348 | XFORM(3, 0, 0, 36); |
349 | |
350 | state->new_tile(state, vs, 4, depth); |
351 | } |
352 | if (depth >= state->max_depth) return 0; |
353 | |
354 | vv_orig = v_trans(v_orig, v_edge); |
355 | |
356 | penrose_p3_large(state, depth+1, -flip, |
357 | vv_orig, v_shrinkphi(v_rotate(v_edge, 180))); |
358 | |
359 | penrose_p3_small(state, depth+1, flip, |
360 | vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip))); |
361 | |
362 | vv_orig = v_trans(v_orig, v_growphi(v_edge)); |
363 | |
364 | penrose_p3_large(state, depth+1, flip, |
365 | vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip))); |
366 | |
367 | |
368 | return 0; |
369 | } |
370 | |
371 | static int penrose_p3_small(penrose_state *state, int depth, int flip, |
372 | vector v_orig, vector v_edge) |
373 | { |
374 | vector vv_orig; |
375 | |
376 | #ifdef DEBUG_PENROSE |
377 | { |
378 | vector vs[3]; |
379 | vs[0] = v_orig; |
380 | XFORM(1, 0, 0, 0); |
381 | XFORM(2, 0, 0, -36*flip); |
382 | |
383 | state->new_tile(state, vs, 3, depth); |
384 | } |
385 | #endif |
386 | |
387 | if (flip > 0) { |
388 | vector vs[4]; |
389 | |
390 | vs[0] = v_orig; |
391 | XFORM(1, 0, 0, -36); |
392 | XFORM(3, 0, 0, 0); |
393 | XFORM(2, 3, 0, -36); |
394 | |
395 | state->new_tile(state, vs, 4, depth); |
396 | } |
397 | if (depth >= state->max_depth) return 0; |
398 | |
399 | /* NB these two are identical to the first two of p3_large. */ |
400 | vv_orig = v_trans(v_orig, v_edge); |
401 | |
402 | penrose_p3_large(state, depth+1, -flip, |
403 | vv_orig, v_shrinkphi(v_rotate(v_edge, 180))); |
404 | |
405 | penrose_p3_small(state, depth+1, flip, |
406 | vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip))); |
407 | |
408 | return 0; |
409 | } |
410 | |
411 | /* ------------------------------------------------------- |
412 | * Utility routines. |
413 | */ |
414 | |
415 | double penrose_side_length(double start_size, int depth) |
416 | { |
417 | return start_size / pow(PHI, depth); |
418 | } |
419 | |
420 | void penrose_count_tiles(int depth, int *nlarge, int *nsmall) |
421 | { |
422 | /* Steal sgt's fibonacci thingummy. */ |
423 | } |
424 | |
425 | /* |
426 | * It turns out that an acute isosceles triangle with sides in ratio 1:phi:phi |
427 | * has an incentre which is conveniently 2*phi^-2 of the way from the apex to |
428 | * the base. Why's that convenient? Because: if we situate the incentre of the |
429 | * triangle at the origin, then we can place the apex at phi^-2 * (B+C), and |
430 | * the other two vertices at apex-B and apex-C respectively. So that's an acute |
431 | * triangle with its long sides of unit length, covering a circle about the |
432 | * origin of radius 1-(2*phi^-2), which is conveniently enough phi^-3. |
433 | * |
434 | * (later mail: this is an overestimate by about 5%) |
435 | */ |
436 | |
a15e5ee3 |
437 | int penrose(penrose_state *state, int which, int angle) |
cebf0b0d |
438 | { |
439 | vector vo = v_origin(); |
440 | vector vb = v_origin(); |
441 | |
442 | vo.b = vo.c = -state->start_size; |
443 | vo = v_shrinkphi(v_shrinkphi(vo)); |
444 | |
445 | vb.b = state->start_size; |
446 | |
a15e5ee3 |
447 | vo = v_rotate(vo, angle); |
448 | vb = v_rotate(vb, angle); |
449 | |
cebf0b0d |
450 | if (which == PENROSE_P2) |
451 | return penrose_p2_large(state, 0, 1, vo, vb); |
452 | else |
453 | return penrose_p3_small(state, 0, 1, vo, vb); |
454 | } |
455 | |
456 | /* |
457 | * We're asked for a MxN grid, which just means a tiling fitting into roughly |
458 | * an MxN space in some kind of reasonable unit - say, the side length of the |
459 | * two-arrow edges of the tiles. By some reasoning in a previous email, that |
460 | * means we want to pick some subarea of a circle of radius 3.11*sqrt(M^2+N^2). |
461 | * To cover that circle, we need to subdivide a triangle large enough that it |
462 | * contains a circle of that radius. |
463 | * |
464 | * Hence: start with those three vectors marking triangle vertices, scale them |
465 | * all up by phi repeatedly until the radius of the inscribed circle gets |
466 | * bigger than the target, and then recurse into that triangle with the same |
467 | * recursion depth as the number of times you scaled up. That will give you |
468 | * tiles of unit side length, covering a circle big enough that if you randomly |
469 | * choose an orientation and coordinates within the circle, you'll be able to |
470 | * get any valid piece of Penrose tiling of size MxN. |
471 | */ |
472 | #define INCIRCLE_RADIUS 0.22426 /* phi^-3 less 5%: see above */ |
473 | |
474 | void penrose_calculate_size(int which, int tilesize, int w, int h, |
475 | double *required_radius, int *start_size, int *depth) |
476 | { |
477 | double rradius, size; |
478 | int n = 0; |
479 | |
480 | /* |
481 | * Fudge factor to scale P2 and P3 tilings differently. This |
482 | * doesn't seem to have much relevance to questions like the |
483 | * average number of tiles per unit area; it's just aesthetic. |
484 | */ |
485 | if (which == PENROSE_P2) |
486 | tilesize = tilesize * 3 / 2; |
487 | else |
488 | tilesize = tilesize * 5 / 4; |
489 | |
490 | rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h)); |
491 | size = tilesize; |
492 | |
493 | while ((size * INCIRCLE_RADIUS) < rradius) { |
494 | n++; |
495 | size = size * PHI; |
496 | } |
497 | |
498 | *start_size = (int)size; |
499 | *depth = n; |
500 | *required_radius = rradius; |
501 | } |
502 | |
503 | /* ------------------------------------------------------- |
504 | * Test code. |
505 | */ |
506 | |
507 | #ifdef TEST_PENROSE |
508 | |
509 | #include <stdio.h> |
510 | #include <string.h> |
511 | |
512 | int show_recursion = 0; |
513 | int ntiles, nfinal; |
514 | |
515 | int test_cb(penrose_state *state, vector *vs, int n, int depth) |
516 | { |
517 | int i, xoff = 0, yoff = 0; |
518 | double l = penrose_side_length(state->start_size, depth); |
519 | double rball = l / 10.0; |
520 | const char *col; |
521 | |
522 | ntiles++; |
523 | if (state->max_depth == depth) { |
524 | col = n == 4 ? "black" : "green"; |
525 | nfinal++; |
526 | } else { |
527 | if (!show_recursion) |
528 | return 0; |
529 | col = n == 4 ? "red" : "blue"; |
530 | } |
531 | if (n != 4) yoff = state->start_size; |
532 | |
533 | printf("<polygon points=\""); |
534 | for (i = 0; i < n; i++) { |
535 | printf("%s%f,%f", (i == 0) ? "" : " ", |
536 | v_x(vs, i) + xoff, v_y(vs, i) + yoff); |
537 | } |
538 | printf("\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n", col, col); |
539 | printf("<ellipse cx=\"%f\" cy=\"%f\" rx=\"%f\" ry=\"%f\" fill=\"%s\" />", |
540 | v_x(vs, 0) + xoff, v_y(vs, 0) + yoff, rball, rball, col); |
541 | |
542 | return 0; |
543 | } |
544 | |
fd66a01d |
545 | void usage_exit(void) |
cebf0b0d |
546 | { |
547 | fprintf(stderr, "Usage: penrose-test [--recursion] P2|P3 SIZE DEPTH\n"); |
548 | exit(1); |
549 | } |
550 | |
551 | int main(int argc, char *argv[]) |
552 | { |
553 | penrose_state ps; |
554 | int which = 0; |
555 | |
556 | while (--argc > 0) { |
557 | char *p = *++argv; |
558 | if (!strcmp(p, "-h") || !strcmp(p, "--help")) { |
559 | usage_exit(); |
560 | } else if (!strcmp(p, "--recursion")) { |
561 | show_recursion = 1; |
562 | } else if (*p == '-') { |
563 | fprintf(stderr, "Unrecognised option '%s'\n", p); |
564 | exit(1); |
565 | } else { |
566 | break; |
567 | } |
568 | } |
569 | |
570 | if (argc < 3) usage_exit(); |
571 | |
572 | if (strcmp(argv[0], "P2") == 0) which = PENROSE_P2; |
573 | else if (strcmp(argv[0], "P3") == 0) which = PENROSE_P3; |
574 | else usage_exit(); |
575 | |
576 | ps.start_size = atoi(argv[1]); |
577 | ps.max_depth = atoi(argv[2]); |
578 | ps.new_tile = test_cb; |
579 | |
580 | ntiles = nfinal = 0; |
581 | |
582 | printf("\ |
583 | <?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n\ |
584 | <!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 20010904//EN\"\n\ |
585 | \"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\ |
586 | \n\ |
587 | <svg xmlns=\"http://www.w3.org/2000/svg\"\n\ |
588 | xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n"); |
589 | |
590 | printf("<g>\n"); |
591 | penrose(&ps, which); |
592 | printf("</g>\n"); |
593 | |
594 | printf("<!-- %d tiles and %d leaf tiles total -->\n", |
595 | ntiles, nfinal); |
596 | |
597 | printf("</svg>"); |
598 | |
599 | return 0; |
600 | } |
601 | |
602 | #endif |
603 | |
604 | #ifdef TEST_VECTORS |
605 | |
606 | static void dbgv(const char *msg, vector v) |
607 | { |
608 | printf("%s: %s\n", msg, v_debug(v)); |
609 | } |
610 | |
611 | int main(int argc, const char *argv[]) |
612 | { |
613 | vector v = v_unit(); |
614 | |
615 | dbgv("unit vector", v); |
616 | v = v_rotate(v, 36); |
617 | dbgv("rotated 36", v); |
618 | v = v_scale(v, 2); |
619 | dbgv("scaled x2", v); |
620 | v = v_shrinkphi(v); |
621 | dbgv("shrunk phi", v); |
622 | v = v_rotate(v, -36); |
623 | dbgv("rotated -36", v); |
624 | |
625 | return 0; |
626 | } |
627 | |
628 | #endif |
629 | /* vim: set shiftwidth=4 tabstop=8: */ |