ref: d35e4b88e85d9bb4f677083c0808b02f51620c2f
dir: /penrose.c/
/* * Generate Penrose tilings via combinatorial coordinates. * * For general explanation of the algorithm: * https://www.chiark.greenend.org.uk/~sgtatham/quasiblog/aperiodic-tilings/ * * I use exactly the same indexing system here that's described in the * article. For the P2 tiling, acute isosceles triangles (half-kites) * are assigned letters A,B, and obtuse ones (half-darts) U,V; for P3, * acute triangles (half of a thin rhomb) are C,D and obtuse ones * (half a thick rhomb) are X,Y. Edges of all triangles are indexed * anticlockwise around the triangle, with 0 being the base and 1,2 * being the two equal legs. */ #include <assert.h> #include <stddef.h> #include <string.h> #include "puzzles.h" #include "penrose.h" #include "penrose-internal.h" #include "tree234.h" bool penrose_valid_letter(char c, int which) { switch (c) { case 'A': case 'B': case 'U': case 'V': return which == PENROSE_P2; case 'C': case 'D': case 'X': case 'Y': return which == PENROSE_P3; default: return false; } } /* * Result of attempting a transition within the coordinate system. * INTERNAL means we've moved to a different child of the same parent, * so the 'internal' substructure gives the type of the new triangle * and which edge of it we came in through; EXTERNAL means we've moved * out of the parent entirely, and the 'external' substructure tells * us which edge of the parent triangle we left by, and if it's * divided in two, which end of that edge (-1 for the left end or +1 * for the right end). If the parent edge is undivided, end == 0. * * The type FAIL _shouldn't_ ever come up! It occurs if you try to * compute an incoming transition with an illegal value of 'end' (i.e. * having the wrong idea of whether the edge is divided), or if you * refer to a child triangle type that doesn't exist in that parent. * If it ever happens in the production code then an assertion will * fail. But it might be useful to other users of the same code. */ typedef struct TransitionResult { enum { INTERNAL, EXTERNAL, FAIL } type; union { struct { char new_child; unsigned char new_edge; } internal; struct { unsigned char parent_edge; signed char end; } external; } u; } TransitionResult; /* Construction function to make an INTERNAL-type TransitionResult */ static inline TransitionResult internal(char new_child, unsigned new_edge) { TransitionResult tr; tr.type = INTERNAL; tr.u.internal.new_child = new_child; tr.u.internal.new_edge = new_edge; return tr; } /* Construction function to make an EXTERNAL-type TransitionResult */ static inline TransitionResult external(unsigned parent_edge, int end) { TransitionResult tr; tr.type = EXTERNAL; tr.u.external.parent_edge = parent_edge; tr.u.external.end = end; return tr; } /* Construction function to make a FAIL-type TransitionResult */ static inline TransitionResult fail(void) { TransitionResult tr; tr.type = FAIL; return tr; } /* * Compute a transition out of a triangle. Can return either INTERNAL * or EXTERNAL types (or FAIL if it gets invalid data). */ static TransitionResult transition(char parent, char child, unsigned edge) { switch (parent) { case 'A': switch (child) { case 'A': switch (edge) { case 0: return external(2, -1); case 1: return external(0, 0); case 2: return internal('B', 1); } case 'B': switch (edge) { case 0: return internal('U', 1); case 1: return internal('A', 2); case 2: return external(1, +1); } case 'U': switch (edge) { case 0: return external(2, +1); case 1: return internal('B', 0); case 2: return external(1, -1); } default: return fail(); } case 'B': switch (child) { case 'A': switch (edge) { case 0: return internal('V', 2); case 1: return external(2, -1); case 2: return internal('B', 1); } case 'B': switch (edge) { case 0: return external(1, +1); case 1: return internal('A', 2); case 2: return external(0, 0); } case 'V': switch (edge) { case 0: return external(1, -1); case 1: return external(2, +1); case 2: return internal('A', 0); } default: return fail(); } case 'U': switch (child) { case 'B': switch (edge) { case 0: return internal('U', 1); case 1: return external(2, 0); case 2: return external(0, +1); } case 'U': switch (edge) { case 0: return external(1, 0); case 1: return internal('B', 0); case 2: return external(0, -1); } default: return fail(); } case 'V': switch (child) { case 'A': switch (edge) { case 0: return internal('V', 2); case 1: return external(0, -1); case 2: return external(1, 0); } case 'V': switch (edge) { case 0: return external(2, 0); case 1: return external(0, +1); case 2: return internal('A', 0); } default: return fail(); } case 'C': switch (child) { case 'C': switch (edge) { case 0: return external(1, +1); case 1: return internal('Y', 1); case 2: return external(0, 0); } case 'Y': switch (edge) { case 0: return external(2, 0); case 1: return internal('C', 1); case 2: return external(1, -1); } default: return fail(); } case 'D': switch (child) { case 'D': switch (edge) { case 0: return external(2, -1); case 1: return external(0, 0); case 2: return internal('X', 2); } case 'X': switch (edge) { case 0: return external(1, 0); case 1: return external(2, +1); case 2: return internal('D', 2); } default: return fail(); } case 'X': switch (child) { case 'C': switch (edge) { case 0: return external(2, +1); case 1: return internal('Y', 1); case 2: return internal('X', 1); } case 'X': switch (edge) { case 0: return external(1, 0); case 1: return internal('C', 2); case 2: return external(0, -1); } case 'Y': switch (edge) { case 0: return external(0, +1); case 1: return internal('C', 1); case 2: return external(2, -1); } default: return fail(); } case 'Y': switch (child) { case 'D': switch (edge) { case 0: return external(1, -1); case 1: return internal('Y', 2); case 2: return internal('X', 2); } case 'X': switch (edge) { case 0: return external(0, -1); case 1: return external(1, +1); case 2: return internal('D', 2); } case 'Y': switch (edge) { case 0: return external(2, 0); case 1: return external(0, +1); case 2: return internal('D', 1); } default: return fail(); } default: return fail(); } } /* * Compute a transition into a parent triangle, after the above * function reported an EXTERNAL transition out of a neighbouring * parent and we had to recurse. Because we're coming inwards, this * should always return an INTERNAL TransitionResult (or FAIL if it * gets invalid data). */ static TransitionResult transition_in(char parent, unsigned edge, int end) { #define EDGEEND(edge, end) (3 * (edge) + 1 + (end)) switch (parent) { case 'A': switch (EDGEEND(edge, end)) { case EDGEEND(0, 0): return internal('A', 1); case EDGEEND(1, -1): return internal('B', 2); case EDGEEND(1, +1): return internal('U', 2); case EDGEEND(2, -1): return internal('U', 0); case EDGEEND(2, +1): return internal('A', 0); default: return fail(); } case 'B': switch (EDGEEND(edge, end)) { case EDGEEND(0, 0): return internal('B', 2); case EDGEEND(1, -1): return internal('B', 0); case EDGEEND(1, +1): return internal('V', 0); case EDGEEND(2, -1): return internal('V', 1); case EDGEEND(2, +1): return internal('A', 1); default: return fail(); } case 'U': switch (EDGEEND(edge, end)) { case EDGEEND(0, -1): return internal('B', 2); case EDGEEND(0, +1): return internal('U', 2); case EDGEEND(1, 0): return internal('U', 0); case EDGEEND(2, 0): return internal('B', 1); default: return fail(); } case 'V': switch (EDGEEND(edge, end)) { case EDGEEND(0, -1): return internal('V', 1); case EDGEEND(0, +1): return internal('A', 1); case EDGEEND(1, 0): return internal('A', 2); case EDGEEND(2, 0): return internal('V', 0); default: return fail(); } case 'C': switch (EDGEEND(edge, end)) { case EDGEEND(0, 0): return internal('C', 2); case EDGEEND(1, -1): return internal('C', 0); case EDGEEND(1, +1): return internal('Y', 2); case EDGEEND(2, 0): return internal('Y', 0); default: return fail(); } case 'D': switch (EDGEEND(edge, end)) { case EDGEEND(0, 0): return internal('D', 1); case EDGEEND(1, 0): return internal('X', 0); case EDGEEND(2, -1): return internal('X', 1); case EDGEEND(2, +1): return internal('D', 0); default: return fail(); } case 'X': switch (EDGEEND(edge, end)) { case EDGEEND(0, -1): return internal('Y', 0); case EDGEEND(0, +1): return internal('X', 2); case EDGEEND(1, 0): return internal('X', 0); case EDGEEND(2, -1): return internal('C', 0); case EDGEEND(2, +1): return internal('Y', 2); default: return fail(); } case 'Y': switch (EDGEEND(edge, end)) { case EDGEEND(0, +1): return internal('X', 0); case EDGEEND(0, -1): return internal('Y', 1); case EDGEEND(1, -1): return internal('X', 1); case EDGEEND(1, +1): return internal('D', 0); case EDGEEND(2, 0): return internal('Y', 0); default: return fail(); } default: return fail(); } #undef EDGEEND } PenroseCoords *penrose_coords_new(void) { PenroseCoords *pc = snew(PenroseCoords); pc->nc = pc->csize = 0; pc->c = NULL; return pc; } void penrose_coords_free(PenroseCoords *pc) { if (pc) { sfree(pc->c); sfree(pc); } } void penrose_coords_make_space(PenroseCoords *pc, size_t size) { if (pc->csize < size) { pc->csize = pc->csize * 5 / 4 + 16; if (pc->csize < size) pc->csize = size; pc->c = sresize(pc->c, pc->csize, char); } } PenroseCoords *penrose_coords_copy(PenroseCoords *pc_in) { PenroseCoords *pc_out = penrose_coords_new(); penrose_coords_make_space(pc_out, pc_in->nc); memcpy(pc_out->c, pc_in->c, pc_in->nc * sizeof(*pc_out->c)); pc_out->nc = pc_in->nc; return pc_out; } /* * The main recursive function for computing the next triangle's * combinatorial coordinates. */ static void penrosectx_step_recurse( PenroseContext *ctx, PenroseCoords *pc, size_t depth, unsigned edge, unsigned *outedge) { TransitionResult tr; penrosectx_extend_coords(ctx, pc, depth+2); /* Look up the transition out of the starting triangle */ tr = transition(pc->c[depth+1], pc->c[depth], edge); /* If we've left the parent triangle, recurse to find out what new * triangle we've landed in at the next size up, and then call * transition_in to find out which child of that parent we're * going to */ if (tr.type == EXTERNAL) { unsigned parent_outedge; penrosectx_step_recurse( ctx, pc, depth+1, tr.u.external.parent_edge, &parent_outedge); tr = transition_in(pc->c[depth+1], parent_outedge, tr.u.external.end); } /* Now we should definitely have ended up in a child of the * (perhaps rewritten) parent triangle */ assert(tr.type == INTERNAL); pc->c[depth] = tr.u.internal.new_child; *outedge = tr.u.internal.new_edge; } void penrosectx_step(PenroseContext *ctx, PenroseCoords *pc, unsigned edge, unsigned *outedge) { /* Allow outedge to be NULL harmlessly, just in case */ unsigned dummy_outedge; if (!outedge) outedge = &dummy_outedge; penrosectx_step_recurse(ctx, pc, 0, edge, outedge); } static Point penrose_triangle_post_edge(char c, unsigned edge) { static const Point acute_post_edge[3] = { {{-1, 1, 0, 1}}, /* phi * t^3 */ {{-1, 1, -1, 1}}, /* t^4 */ {{-1, 1, 0, 0}}, /* 1/phi * t^3 */ }; static const Point obtuse_post_edge[3] = { {{0, -1, 1, 0}}, /* 1/phi * t^4 */ {{0, 0, 1, 0}}, /* t^2 */ {{-1, 0, 0, 1}}, /* phi * t^4 */ }; switch (c) { case 'A': case 'B': case 'C': case 'D': return acute_post_edge[edge]; default: /* case 'U': case 'V': case 'X': case 'Y': */ return obtuse_post_edge[edge]; } } void penrose_place(PenroseTriangle *tri, Point u, Point v, int index_of_u) { unsigned i; Point here = u, delta = point_sub(v, u); for (i = 0; i < 3; i++) { unsigned edge = (index_of_u + i) % 3; tri->vertices[edge] = here; here = point_add(here, delta); delta = point_mul(delta, penrose_triangle_post_edge( tri->pc->c[0], edge)); } } void penrose_free(PenroseTriangle *tri) { penrose_coords_free(tri->pc); sfree(tri); } static bool penrose_relative_probability(char c) { /* Penrose tile probability ratios are always phi, so we can * approximate that very well using two consecutive Fibonacci * numbers */ switch (c) { case 'A': case 'B': case 'X': case 'Y': return 165580141; case 'C': case 'D': case 'U': case 'V': return 102334155; default: return 0; } } static char penrose_choose_random(const char *possibilities, random_state *rs) { const char *p; unsigned long value, limit = 0; for (p = possibilities; *p; p++) limit += penrose_relative_probability(*p); value = random_upto(rs, limit); for (p = possibilities; *p; p++) { unsigned long curr = penrose_relative_probability(*p); if (value < curr) return *p; value -= curr; } assert(false && "Probability overflow!"); return possibilities[0]; } static const char *penrose_starting_tiles(int which) { return which == PENROSE_P2 ? "ABUV" : "CDXY"; } static const char *penrose_valid_parents(char tile) { switch (tile) { case 'A': return "ABV"; case 'B': return "ABU"; case 'U': return "AU"; case 'V': return "BV"; case 'C': return "CX"; case 'D': return "DY"; case 'X': return "DXY"; case 'Y': return "CXY"; default: return NULL; } } void penrosectx_init_random(PenroseContext *ctx, random_state *rs, int which) { ctx->rs = rs; ctx->must_free_rs = false; ctx->prototype = penrose_coords_new(); penrose_coords_make_space(ctx->prototype, 1); ctx->prototype->c[0] = penrose_choose_random( penrose_starting_tiles(which), rs); ctx->prototype->nc = 1; ctx->start_vertex = random_upto(rs, 3); ctx->orientation = random_upto(rs, 10); } void penrosectx_init_from_params( PenroseContext *ctx, const struct PenrosePatchParams *ps) { size_t i; ctx->rs = NULL; ctx->must_free_rs = false; ctx->prototype = penrose_coords_new(); penrose_coords_make_space(ctx->prototype, ps->ncoords); for (i = 0; i < ps->ncoords; i++) ctx->prototype->c[i] = ps->coords[i]; ctx->prototype->nc = ps->ncoords; ctx->start_vertex = ps->start_vertex; ctx->orientation = ps->orientation; } void penrosectx_cleanup(PenroseContext *ctx) { if (ctx->must_free_rs) random_free(ctx->rs); penrose_coords_free(ctx->prototype); } PenroseCoords *penrosectx_initial_coords(PenroseContext *ctx) { return penrose_coords_copy(ctx->prototype); } void penrosectx_extend_coords(PenroseContext *ctx, PenroseCoords *pc, size_t n) { if (ctx->prototype->nc < n) { penrose_coords_make_space(ctx->prototype, n); while (ctx->prototype->nc < n) { if (!ctx->rs) { /* * For safety, similarly to spectre.c, we respond to a * lack of available random_state by making a * deterministic one. */ ctx->rs = random_new("dummy", 5); ctx->must_free_rs = true; } ctx->prototype->c[ctx->prototype->nc] = penrose_choose_random( penrose_valid_parents(ctx->prototype->c[ctx->prototype->nc-1]), ctx->rs); ctx->prototype->nc++; } } penrose_coords_make_space(pc, n); while (pc->nc < n) { pc->c[pc->nc] = ctx->prototype->c[pc->nc]; pc->nc++; } } static Point penrose_triangle_edge_0_length(char c) { static const Point one = {{ 1, 0, 0, 0 }}; static const Point phi = {{ 1, 0, 1, -1 }}; static const Point invphi = {{ 0, 0, 1, -1 }}; switch (c) { /* P2 tiling: unit-length edges are the long edges, i.e. edges * 1,2 of AB and edge 0 of UV. So AB have edge 0 short. */ case 'A': case 'B': return invphi; case 'U': case 'V': return one; /* P3 tiling: unit-length edges are edges 1,2 of everything, * so CD have edge 0 short and XY have it long. */ case 'C': case 'D': return invphi; default: /* case 'X': case 'Y': */ return phi; } } PenroseTriangle *penrose_initial(PenroseContext *ctx) { char type = ctx->prototype->c[0]; Point origin = {{ 0, 0, 0, 0 }}; Point edge0 = penrose_triangle_edge_0_length(type); Point negoffset; size_t i; PenroseTriangle *tri = snew(PenroseTriangle); /* Orient the triangle by deciding what vector edge #0 should traverse */ edge0 = point_mul(edge0, point_rot(ctx->orientation)); /* Place the triangle at an arbitrary position, in that orientation */ tri->pc = penrose_coords_copy(ctx->prototype); penrose_place(tri, origin, edge0, 0); /* Now translate so that the appropriate vertex is at the origin */ negoffset = tri->vertices[ctx->start_vertex]; for (i = 0; i < 3; i++) tri->vertices[i] = point_sub(tri->vertices[i], negoffset); return tri; } PenroseTriangle *penrose_adjacent(PenroseContext *ctx, const PenroseTriangle *src_tri, unsigned src_edge, unsigned *dst_edge_out) { unsigned dst_edge; PenroseTriangle *dst_tri = snew(PenroseTriangle); dst_tri->pc = penrose_coords_copy(src_tri->pc); penrosectx_step(ctx, dst_tri->pc, src_edge, &dst_edge); penrose_place(dst_tri, src_tri->vertices[(src_edge+1) % 3], src_tri->vertices[src_edge], dst_edge); if (dst_edge_out) *dst_edge_out = dst_edge; return dst_tri; } static int penrose_cmp(void *av, void *bv) { PenroseTriangle *a = (PenroseTriangle *)av, *b = (PenroseTriangle *)bv; size_t i, j; /* We should only ever need to compare the first two vertices of * any triangle, because those force the rest */ for (i = 0; i < 2; i++) { for (j = 0; j < 4; j++) { int ac = a->vertices[i].coeffs[j], bc = b->vertices[i].coeffs[j]; if (ac < bc) return -1; if (ac > bc) return +1; } } return 0; } static unsigned penrose_sibling_edge_index(char c) { switch (c) { case 'A': case 'U': return 2; case 'B': case 'V': return 1; default: return 0; } } void penrosectx_generate( PenroseContext *ctx, bool (*inbounds)(void *inboundsctx, const PenroseTriangle *tri), void *inboundsctx, void (*tile)(void *tilectx, const Point *vertices), void *tilectx) { tree234 *placed = newtree234(penrose_cmp); PenroseTriangle *qhead = NULL, *qtail = NULL; { PenroseTriangle *tri = penrose_initial(ctx); add234(placed, tri); tri->next = NULL; tri->reported = false; if (inbounds(inboundsctx, tri)) qhead = qtail = tri; } while (qhead) { PenroseTriangle *tri = qhead; unsigned edge; unsigned sibling_edge = penrose_sibling_edge_index(tri->pc->c[0]); for (edge = 0; edge < 3; edge++) { PenroseTriangle *new_tri, *found_tri; new_tri = penrose_adjacent(ctx, tri, edge, NULL); if (!inbounds(inboundsctx, new_tri)) { penrose_free(new_tri); continue; } found_tri = find234(placed, new_tri, NULL); if (found_tri) { if (edge == sibling_edge && !tri->reported && !found_tri->reported) { /* * found_tri and tri are opposite halves of the * same tile; both are in the tree, and haven't * yet been reported as a completed tile. */ unsigned new_sibling_edge = penrose_sibling_edge_index( found_tri->pc->c[0]); Point tilevertices[4] = { tri->vertices[(sibling_edge + 1) % 3], tri->vertices[(sibling_edge + 2) % 3], found_tri->vertices[(new_sibling_edge + 1) % 3], found_tri->vertices[(new_sibling_edge + 2) % 3], }; tile(tilectx, tilevertices); tri->reported = true; found_tri->reported = true; } penrose_free(new_tri); continue; } add234(placed, new_tri); qtail->next = new_tri; qtail = new_tri; new_tri->next = NULL; new_tri->reported = false; } qhead = qhead->next; } { PenroseTriangle *tri; while ((tri = delpos234(placed, 0)) != NULL) penrose_free(tri); freetree234(placed); } } const char *penrose_tiling_params_invalid( const struct PenrosePatchParams *params, int which) { size_t i; if (params->ncoords == 0) return "expected at least one coordinate"; for (i = 0; i < params->ncoords; i++) { if (!penrose_valid_letter(params->coords[i], which)) return "invalid coordinate letter"; if (i > 0 && !strchr(penrose_valid_parents(params->coords[i-1]), params->coords[i])) return "invalid pair of consecutive coordinates"; } return NULL; } struct PenroseOutputCtx { int xoff, yoff; Coord xmin, xmax, ymin, ymax; penrose_tile_callback_fn external_cb; void *external_cbctx; }; static bool inbounds(void *vctx, const PenroseTriangle *tri) { struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx; size_t i; for (i = 0; i < 3; i++) { Coord x = point_x(tri->vertices[i]); Coord y = point_y(tri->vertices[i]); if (coord_cmp(x, octx->xmin) < 0 || coord_cmp(x, octx->xmax) > 0 || coord_cmp(y, octx->ymin) < 0 || coord_cmp(y, octx->ymax) > 0) return false; } return true; } static void null_output_tile(void *vctx, const Point *vertices) { } static void really_output_tile(void *vctx, const Point *vertices) { struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx; size_t i; int coords[16]; for (i = 0; i < 4; i++) { Coord x = point_x(vertices[i]); Coord y = point_y(vertices[i]); coords[4*i + 0] = x.c1 + octx->xoff; coords[4*i + 1] = x.cr5; coords[4*i + 2] = y.c1 + octx->yoff; coords[4*i + 3] = y.cr5; } octx->external_cb(octx->external_cbctx, coords); } static void penrose_set_bounds(struct PenroseOutputCtx *octx, int w, int h) { octx->xoff = w/2; octx->yoff = h/2; octx->xmin.c1 = -octx->xoff; octx->xmax.c1 = -octx->xoff + w; octx->ymin.c1 = octx->yoff - h; octx->ymax.c1 = octx->yoff; octx->xmin.cr5 = 0; octx->xmax.cr5 = 0; octx->ymin.cr5 = 0; octx->ymax.cr5 = 0; } void penrose_tiling_randomise(struct PenrosePatchParams *params, int which, int w, int h, random_state *rs) { PenroseContext ctx[1]; struct PenroseOutputCtx octx[1]; penrose_set_bounds(octx, w, h); penrosectx_init_random(ctx, rs, which); penrosectx_generate(ctx, inbounds, octx, null_output_tile, NULL); params->orientation = ctx->orientation; params->start_vertex = ctx->start_vertex; params->ncoords = ctx->prototype->nc; params->coords = snewn(params->ncoords, char); memcpy(params->coords, ctx->prototype->c, params->ncoords); penrosectx_cleanup(ctx); } void penrose_tiling_generate( const struct PenrosePatchParams *params, int w, int h, penrose_tile_callback_fn cb, void *cbctx) { PenroseContext ctx[1]; struct PenroseOutputCtx octx[1]; penrose_set_bounds(octx, w, h); octx->external_cb = cb; octx->external_cbctx = cbctx; penrosectx_init_from_params(ctx, params); penrosectx_generate(ctx, inbounds, octx, really_output_tile, octx); penrosectx_cleanup(ctx); }