chiark / gitweb /
Tents: mark squares as non-tents with {Shift,Control}-cursor keys.
[sgt-puzzles.git] / grid.c
diff --git a/grid.c b/grid.c
index 1415e61f067a7b520b40f5f30826bef278400c25..6a90c9942e3a44edff78696e9c6ffd00e554cdce 100644 (file)
--- a/grid.c
+++ b/grid.c
 #include <assert.h>
 #include <ctype.h>
 #include <math.h>
+#include <float.h>
 
 #include "puzzles.h"
 #include "tree234.h"
 #include "grid.h"
+#include "penrose.h"
 
 /* Debugging options */
 
@@ -50,14 +52,13 @@ void grid_free(grid *g)
 
 /* Used by the other grid generators.  Create a brand new grid with nothing
  * initialised (all lists are NULL) */
-static grid *grid_new(void)
+static grid *grid_empty(void)
 {
     grid *g = snew(grid);
     g->faces = NULL;
     g->edges = NULL;
     g->dots = NULL;
     g->num_faces = g->num_edges = g->num_dots = 0;
-    g->middle_face = NULL;
     g->refcount = 1;
     g->lowest_x = g->lowest_y = g->highest_x = g->highest_y = 0;
     return g;
@@ -92,17 +93,9 @@ static double point_line_distance(long px, long py,
  * Returns the nearest edge, or NULL if no edge is reasonably
  * near the position.
  *
- * This algorithm is nice and generic, and doesn't depend on any particular
- * geometric layout of the grid:
- *   Start at any dot (pick one next to middle_face).
- *   Walk along a path by choosing, from all nearby dots, the one that is
- *   nearest the target (x,y).  Hopefully end up at the dot which is closest
- *   to (x,y).  Should work, as long as faces aren't too badly shaped.
- *   Then examine each edge around this dot, and pick whichever one is
- *   closest (perpendicular distance) to (x,y).
- *   Using perpendicular distance is not quite right - the edge might be
- *   "off to one side".  So we insist that the triangle with (x,y) has
- *   acute angles at the edge's dots.
+ * Just judging edges by perpendicular distance is not quite right -
+ * the edge might be "off to one side". So we insist that the triangle
+ * with (x,y) has acute angles at the edge's dots.
  *
  *     edge1
  *  *---------*------
@@ -116,50 +109,14 @@ static double point_line_distance(long px, long py,
  */
 grid_edge *grid_nearest_edge(grid *g, int x, int y)
 {
-    grid_dot *cur;
     grid_edge *best_edge;
     double best_distance = 0;
     int i;
 
-    cur = g->middle_face->dots[0];
-
-    for (;;) {
-        /* Target to beat */
-        long dist = SQ((long)cur->x - (long)x) + SQ((long)cur->y - (long)y);
-        /* Look for nearer dot - if found, store in 'new'. */
-        grid_dot *new = cur;
-        int i;
-        /* Search all dots in all faces touching this dot.  Some shapes
-         * (such as in Cairo) don't quite work properly if we only search
-         * the dot's immediate neighbours. */
-        for (i = 0; i < cur->order; i++) {
-            grid_face *f = cur->faces[i];
-            int j;
-            if (!f) continue;
-            for (j = 0; j < f->order; j++) {
-               long new_dist;
-                grid_dot *d = f->dots[j];
-                if (d == cur) continue;
-                new_dist = SQ((long)d->x - (long)x) + SQ((long)d->y - (long)y);
-                if (new_dist < dist) { /* found closer dot */
-                    new = d;
-                    dist = new_dist;
-                }
-            }
-        }
-
-        if (new == cur) {
-            /* Didn't find a closer dot among the neighbours of 'cur' */
-            break;
-        } else {
-            cur = new;
-        }
-    }
-    /* 'cur' is nearest dot, so find which of the dot's edges is closest. */
     best_edge = NULL;
 
-    for (i = 0; i < cur->order; i++) {
-        grid_edge *e = cur->edges[i];
+    for (i = 0; i < g->num_edges; i++) {
+        grid_edge *e = &g->edges[i];
         long e2; /* squared length of edge */
         long a2, b2; /* squared lengths of other sides */
         double dist;
@@ -203,13 +160,95 @@ grid_edge *grid_nearest_edge(grid *g, int x, int y)
  * Grid generation
  */
 
-#ifdef DEBUG_GRID
+#ifdef SVG_GRID
+
+#define SVG_DOTS  1
+#define SVG_EDGES 2
+#define SVG_FACES 4
+
+#define FACE_COLOUR "red"
+#define EDGE_COLOUR "blue"
+#define DOT_COLOUR "black"
+
+static void grid_output_svg(FILE *fp, grid *g, int which)
+{
+    int i, j;
+
+    fprintf(fp,"\
+<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n\
+<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 20010904//EN\"\n\
+\"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\
+\n\
+<svg xmlns=\"http://www.w3.org/2000/svg\"\n\
+xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n");
+
+    if (which & SVG_FACES) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_faces; i++) {
+            grid_face *f = g->faces + i;
+            fprintf(fp, "<polygon points=\"");
+            for (j = 0; j < f->order; j++) {
+                grid_dot *d = f->dots[j];
+                fprintf(fp, "%s%d,%d", (j == 0) ? "" : " ",
+                        d->x, d->y);
+            }
+            fprintf(fp, "\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n",
+                    FACE_COLOUR, FACE_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+    if (which & SVG_EDGES) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_edges; i++) {
+            grid_edge *e = g->edges + i;
+            grid_dot *d1 = e->dot1, *d2 = e->dot2;
+
+            fprintf(fp, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" "
+                        "style=\"stroke: %s\" />\n",
+                        d1->x, d1->y, d2->x, d2->y, EDGE_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+
+    if (which & SVG_DOTS) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_dots; i++) {
+            grid_dot *d = g->dots + i;
+            fprintf(fp, "<ellipse cx=\"%d\" cy=\"%d\" rx=\"%d\" ry=\"%d\" fill=\"%s\" />",
+                    d->x, d->y, g->tilesize/20, g->tilesize/20, DOT_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+
+    fprintf(fp, "</svg>\n");
+}
+#endif
+
+#ifdef SVG_GRID
+#include <errno.h>
+
+static void grid_try_svg(grid *g, int which)
+{
+    char *svg = getenv("PUZZLES_SVG_GRID");
+    if (svg) {
+        FILE *svgf = fopen(svg, "w");
+        if (svgf) {
+            grid_output_svg(svgf, g, which);
+            fclose(svgf);
+        } else {
+            fprintf(stderr, "Unable to open file `%s': %s", svg, strerror(errno));
+        }
+    }
+}
+#endif
+
 /* Show the basic grid information, before doing grid_make_consistent */
-static void grid_print_basic(grid *g)
+static void grid_debug_basic(grid *g)
 {
     /* TODO: Maybe we should generate an SVG image of the dots and lines
      * of the grid here, before grid_make_consistent.
      * Would help with debugging grid generation. */
+#ifdef DEBUG_GRID
     int i;
     printf("--- Basic Grid Data ---\n");
     for (i = 0; i < g->num_faces; i++) {
@@ -222,11 +261,16 @@ static void grid_print_basic(grid *g)
         }
         printf("]\n");
     }
-    printf("Middle face: %d\n", (int)(g->middle_face - g->faces));
+#endif
+#ifdef SVG_GRID
+    grid_try_svg(g, SVG_FACES);
+#endif
 }
+
 /* Show the derived grid information, computed by grid_make_consistent */
-static void grid_print_derived(grid *g)
+static void grid_debug_derived(grid *g)
 {
+#ifdef DEBUG_GRID
     /* edges */
     int i;
     printf("--- Derived Grid Data ---\n");
@@ -266,8 +310,11 @@ static void grid_print_derived(grid *g)
         }
         printf("]\n");
     }
+#endif
+#ifdef SVG_GRID
+    grid_try_svg(g, SVG_DOTS | SVG_EDGES | SVG_FACES);
+#endif
 }
-#endif /* DEBUG_GRID */
 
 /* Helper function for building incomplete-edges list in
  * grid_make_consistent() */
@@ -296,6 +343,155 @@ static int grid_edge_bydots_cmpfn(void *v1, void *v2)
     return 0;
 }
 
+/*
+ * 'Vigorously trim' a grid, by which I mean deleting any isolated or
+ * uninteresting faces. By which, in turn, I mean: ensure that the
+ * grid is composed solely of faces adjacent to at least one
+ * 'landlocked' dot (i.e. one not in contact with the infinite
+ * exterior face), and that all those dots are in a single connected
+ * component.
+ *
+ * This function operates on, and returns, a grid satisfying the
+ * preconditions to grid_make_consistent() rather than the
+ * postconditions. (So call it first.)
+ */
+static void grid_trim_vigorously(grid *g)
+{
+    int *dotpairs, *faces, *dots;
+    int *dsf;
+    int i, j, k, size, newfaces, newdots;
+
+    /*
+     * First construct a matrix in which each ordered pair of dots is
+     * mapped to the index of the face in which those dots occur in
+     * that order.
+     */
+    dotpairs = snewn(g->num_dots * g->num_dots, int);
+    for (i = 0; i < g->num_dots; i++)
+        for (j = 0; j < g->num_dots; j++)
+            dotpairs[i*g->num_dots+j] = -1;
+    for (i = 0; i < g->num_faces; i++) {
+        grid_face *f = g->faces + i;
+        int dot0 = f->dots[f->order-1] - g->dots;
+        for (j = 0; j < f->order; j++) {
+            int dot1 = f->dots[j] - g->dots;
+            dotpairs[dot0 * g->num_dots + dot1] = i;
+            dot0 = dot1;
+        }
+    }
+
+    /*
+     * Now we can identify landlocked dots: they're the ones all of
+     * whose edges have a mirror-image counterpart in this matrix.
+     */
+    dots = snewn(g->num_dots, int);
+    for (i = 0; i < g->num_dots; i++) {
+        dots[i] = TRUE;
+        for (j = 0; j < g->num_dots; j++) {
+            if ((dotpairs[i*g->num_dots+j] >= 0) ^
+                (dotpairs[j*g->num_dots+i] >= 0))
+                dots[i] = FALSE;    /* non-duplicated edge: coastal dot */
+        }
+    }
+
+    /*
+     * Now identify connected pairs of landlocked dots, and form a dsf
+     * unifying them.
+     */
+    dsf = snew_dsf(g->num_dots);
+    for (i = 0; i < g->num_dots; i++)
+        for (j = 0; j < i; j++)
+            if (dots[i] && dots[j] &&
+                dotpairs[i*g->num_dots+j] >= 0 &&
+                dotpairs[j*g->num_dots+i] >= 0)
+                dsf_merge(dsf, i, j);
+
+    /*
+     * Now look for the largest component.
+     */
+    size = 0;
+    j = -1;
+    for (i = 0; i < g->num_dots; i++) {
+        int newsize;
+        if (dots[i] && dsf_canonify(dsf, i) == i &&
+            (newsize = dsf_size(dsf, i)) > size) {
+            j = i;
+            size = newsize;
+        }
+    }
+
+    /*
+     * Work out which faces we're going to keep (precisely those with
+     * at least one dot in the same connected component as j) and
+     * which dots (those required by any face we're keeping).
+     *
+     * At this point we reuse the 'dots' array to indicate the dots
+     * we're keeping, rather than the ones that are landlocked.
+     */
+    faces = snewn(g->num_faces, int);
+    for (i = 0; i < g->num_faces; i++)
+        faces[i] = 0;
+    for (i = 0; i < g->num_dots; i++)
+        dots[i] = 0;
+    for (i = 0; i < g->num_faces; i++) {
+        grid_face *f = g->faces + i;
+        int keep = FALSE;
+        for (k = 0; k < f->order; k++)
+            if (dsf_canonify(dsf, f->dots[k] - g->dots) == j)
+                keep = TRUE;
+        if (keep) {
+            faces[i] = TRUE;
+            for (k = 0; k < f->order; k++)
+                dots[f->dots[k]-g->dots] = TRUE;
+        }
+    }
+
+    /*
+     * Work out the new indices of those faces and dots, when we
+     * compact the arrays containing them.
+     */
+    for (i = newfaces = 0; i < g->num_faces; i++)
+        faces[i] = (faces[i] ? newfaces++ : -1);
+    for (i = newdots = 0; i < g->num_dots; i++)
+        dots[i] = (dots[i] ? newdots++ : -1);
+
+    /*
+     * Free the dynamically allocated 'dots' pointer lists in faces
+     * we're going to discard.
+     */
+    for (i = 0; i < g->num_faces; i++)
+        if (faces[i] < 0)
+            sfree(g->faces[i].dots);
+
+    /*
+     * Go through and compact the arrays.
+     */
+    for (i = 0; i < g->num_dots; i++)
+        if (dots[i] >= 0) {
+            grid_dot *dnew = g->dots + dots[i], *dold = g->dots + i;
+            *dnew = *dold;             /* structure copy */
+        }
+    for (i = 0; i < g->num_faces; i++)
+        if (faces[i] >= 0) {
+            grid_face *fnew = g->faces + faces[i], *fold = g->faces + i;
+            *fnew = *fold;             /* structure copy */
+            for (j = 0; j < fnew->order; j++) {
+                /*
+                 * Reindex the dots in this face.
+                 */
+                k = fnew->dots[j] - g->dots;
+                fnew->dots[j] = g->dots + dots[k];
+            }
+        }
+    g->num_faces = newfaces;
+    g->num_dots = newdots;
+
+    sfree(dotpairs);
+    sfree(dsf);
+    sfree(dots);
+    sfree(faces);
+}
+
 /* Input: grid has its dots and faces initialised:
  * - dots have (optionally) x and y coordinates, but no edges or faces
  * (pointers are NULL).
@@ -311,9 +507,7 @@ static void grid_make_consistent(grid *g)
     tree234 *incomplete_edges;
     grid_edge *next_new_edge; /* Where new edge will go into g->edges */
 
-#ifdef DEBUG_GRID
-    grid_print_basic(g);
-#endif
+    grid_debug_basic(g);
 
     /* ====== Stage 1 ======
      * Generate edges
@@ -591,10 +785,8 @@ static void grid_make_consistent(grid *g)
             g->highest_y = max(g->highest_y, d->y);
         }
     }
-    
-#ifdef DEBUG_GRID
-    grid_print_derived(g);
-#endif
+
+    grid_debug_derived(g);
 }
 
 /* Helpers for making grid-generation easier.  These functions are only
@@ -621,6 +813,7 @@ static void grid_face_add_new(grid *g, int face_size)
     for (i = 0; i < face_size; i++)
         new_face->dots[i] = NULL;
     new_face->edges = NULL;
+    new_face->has_incentre = FALSE;
     g->num_faces++;
 }
 /* Assumes dot list has enough space */
@@ -666,6 +859,516 @@ static void grid_face_set_dot(grid *g, grid_dot *d, int position)
     last_face->dots[position] = d;
 }
 
+/*
+ * Helper routines for grid_find_incentre.
+ */
+static int solve_2x2_matrix(double mx[4], double vin[2], double vout[2])
+{
+    double inv[4];
+    double det;
+    det = (mx[0]*mx[3] - mx[1]*mx[2]);
+    if (det == 0)
+        return FALSE;
+
+    inv[0] = mx[3] / det;
+    inv[1] = -mx[1] / det;
+    inv[2] = -mx[2] / det;
+    inv[3] = mx[0] / det;
+
+    vout[0] = inv[0]*vin[0] + inv[1]*vin[1];
+    vout[1] = inv[2]*vin[0] + inv[3]*vin[1];
+
+    return TRUE;
+}
+static int solve_3x3_matrix(double mx[9], double vin[3], double vout[3])
+{
+    double inv[9];
+    double det;
+
+    det = (mx[0]*mx[4]*mx[8] + mx[1]*mx[5]*mx[6] + mx[2]*mx[3]*mx[7] -
+           mx[0]*mx[5]*mx[7] - mx[1]*mx[3]*mx[8] - mx[2]*mx[4]*mx[6]);
+    if (det == 0)
+        return FALSE;
+
+    inv[0] = (mx[4]*mx[8] - mx[5]*mx[7]) / det;
+    inv[1] = (mx[2]*mx[7] - mx[1]*mx[8]) / det;
+    inv[2] = (mx[1]*mx[5] - mx[2]*mx[4]) / det;
+    inv[3] = (mx[5]*mx[6] - mx[3]*mx[8]) / det;
+    inv[4] = (mx[0]*mx[8] - mx[2]*mx[6]) / det;
+    inv[5] = (mx[2]*mx[3] - mx[0]*mx[5]) / det;
+    inv[6] = (mx[3]*mx[7] - mx[4]*mx[6]) / det;
+    inv[7] = (mx[1]*mx[6] - mx[0]*mx[7]) / det;
+    inv[8] = (mx[0]*mx[4] - mx[1]*mx[3]) / det;
+
+    vout[0] = inv[0]*vin[0] + inv[1]*vin[1] + inv[2]*vin[2];
+    vout[1] = inv[3]*vin[0] + inv[4]*vin[1] + inv[5]*vin[2];
+    vout[2] = inv[6]*vin[0] + inv[7]*vin[1] + inv[8]*vin[2];
+
+    return TRUE;
+}
+
+void grid_find_incentre(grid_face *f)
+{
+    double xbest, ybest, bestdist;
+    int i, j, k, m;
+    grid_dot *edgedot1[3], *edgedot2[3];
+    grid_dot *dots[3];
+    int nedges, ndots;
+
+    if (f->has_incentre)
+        return;
+
+    /*
+     * Find the point in the polygon with the maximum distance to any
+     * edge or corner.
+     *
+     * Such a point must exist which is in contact with at least three
+     * edges and/or vertices. (Proof: if it's only in contact with two
+     * edges and/or vertices, it can't even be at a _local_ maximum -
+     * any such circle can always be expanded in some direction.) So
+     * we iterate through all 3-subsets of the combined set of edges
+     * and vertices; for each subset we generate one or two candidate
+     * points that might be the incentre, and then we vet each one to
+     * see if it's inside the polygon and what its maximum radius is.
+     *
+     * (There's one case which this algorithm will get noticeably
+     * wrong, and that's when a continuum of equally good answers
+     * exists due to parallel edges. Consider a long thin rectangle,
+     * for instance, or a parallelogram. This algorithm will pick a
+     * point near one end, and choose the end arbitrarily; obviously a
+     * nicer point to choose would be in the centre. To fix this I
+     * would have to introduce a special-case system which detected
+     * parallel edges in advance, set aside all candidate points
+     * generated using both edges in a parallel pair, and generated
+     * some additional candidate points half way between them. Also,
+     * of course, I'd have to cope with rounding error making such a
+     * point look worse than one of its endpoints. So I haven't done
+     * this for the moment, and will cross it if necessary when I come
+     * to it.)
+     *
+     * We don't actually iterate literally over _edges_, in the sense
+     * of grid_edge structures. Instead, we fill in edgedot1[] and
+     * edgedot2[] with a pair of dots adjacent in the face's list of
+     * vertices. This ensures that we get the edges in consistent
+     * orientation, which we could not do from the grid structure
+     * alone. (A moment's consideration of an order-3 vertex should
+     * make it clear that if a notional arrow was written on each
+     * edge, _at least one_ of the three faces bordering that vertex
+     * would have to have the two arrows tip-to-tip or tail-to-tail
+     * rather than tip-to-tail.)
+     */
+    nedges = ndots = 0;
+    bestdist = 0;
+    xbest = ybest = 0;
+
+    for (i = 0; i+2 < 2*f->order; i++) {
+        if (i < f->order) {
+            edgedot1[nedges] = f->dots[i];
+            edgedot2[nedges++] = f->dots[(i+1)%f->order];
+        } else
+            dots[ndots++] = f->dots[i - f->order];
+
+        for (j = i+1; j+1 < 2*f->order; j++) {
+            if (j < f->order) {
+                edgedot1[nedges] = f->dots[j];
+                edgedot2[nedges++] = f->dots[(j+1)%f->order];
+            } else
+                dots[ndots++] = f->dots[j - f->order];
+
+            for (k = j+1; k < 2*f->order; k++) {
+                double cx[2], cy[2];   /* candidate positions */
+                int cn = 0;            /* number of candidates */
+
+                if (k < f->order) {
+                    edgedot1[nedges] = f->dots[k];
+                    edgedot2[nedges++] = f->dots[(k+1)%f->order];
+                } else
+                    dots[ndots++] = f->dots[k - f->order];
+
+                /*
+                 * Find a point, or pair of points, equidistant from
+                 * all the specified edges and/or vertices.
+                 */
+                if (nedges == 3) {
+                    /*
+                     * Three edges. This is a linear matrix equation:
+                     * each row of the matrix represents the fact that
+                     * the point (x,y) we seek is at distance r from
+                     * that edge, and we solve three of those
+                     * simultaneously to obtain x,y,r. (We ignore r.)
+                     */
+                    double matrix[9], vector[3], vector2[3];
+                    int m;
+
+                    for (m = 0; m < 3; m++) {
+                        int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
+                        int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
+                        int dx = x2-x1, dy = y2-y1;
+
+                        /*
+                         * ((x,y) - (x1,y1)) . (dy,-dx) = r |(dx,dy)|
+                         *
+                         * => x dy - y dx - r |(dx,dy)| = (x1 dy - y1 dx)
+                         */
+                        matrix[3*m+0] = dy;
+                        matrix[3*m+1] = -dx;
+                        matrix[3*m+2] = -sqrt((double)dx*dx+(double)dy*dy);
+                        vector[m] = (double)x1*dy - (double)y1*dx;
+                    }
+
+                    if (solve_3x3_matrix(matrix, vector, vector2)) {
+                        cx[cn] = vector2[0];
+                        cy[cn] = vector2[1];
+                        cn++;
+                    }
+                } else if (nedges == 2) {
+                    /*
+                     * Two edges and a dot. This will end up in a
+                     * quadratic equation.
+                     *
+                     * First, look at the two edges. Having our point
+                     * be some distance r from both of them gives rise
+                     * to a pair of linear equations in x,y,r of the
+                     * form
+                     *
+                     *   (x-x1) dy - (y-y1) dx = r sqrt(dx^2+dy^2)
+                     *
+                     * We eliminate r between those equations to give
+                     * us a single linear equation in x,y describing
+                     * the locus of points equidistant from both lines
+                     * - i.e. the angle bisector. 
+                     *
+                     * We then choose one of x,y to be a parameter t,
+                     * and derive linear formulae for x,y,r in terms
+                     * of t. This enables us to write down the
+                     * circular equation (x-xd)^2+(y-yd)^2=r^2 as a
+                     * quadratic in t; solving that and substituting
+                     * in for x,y gives us two candidate points.
+                     */
+                    double eqs[2][4];  /* a,b,c,d : ax+by+cr=d */
+                    double eq[3];      /* a,b,c: ax+by=c */
+                    double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
+                    double q[3];                /* a,b,c: at^2+bt+c=0 */
+                    double disc;
+
+                    /* Find equations of the two input lines. */
+                    for (m = 0; m < 2; m++) {
+                        int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
+                        int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
+                        int dx = x2-x1, dy = y2-y1;
+
+                        eqs[m][0] = dy;
+                        eqs[m][1] = -dx;
+                        eqs[m][2] = -sqrt(dx*dx+dy*dy);
+                        eqs[m][3] = x1*dy - y1*dx;
+                    }
+
+                    /* Derive the angle bisector by eliminating r. */
+                    eq[0] = eqs[0][0]*eqs[1][2] - eqs[1][0]*eqs[0][2];
+                    eq[1] = eqs[0][1]*eqs[1][2] - eqs[1][1]*eqs[0][2];
+                    eq[2] = eqs[0][3]*eqs[1][2] - eqs[1][3]*eqs[0][2];
+
+                    /* Parametrise x and y in terms of some t. */
+                    if (fabs(eq[0]) < fabs(eq[1])) {
+                        /* Parameter is x. */
+                        xt[0] = 1; xt[1] = 0;
+                        yt[0] = -eq[0]/eq[1]; yt[1] = eq[2]/eq[1];
+                    } else {
+                        /* Parameter is y. */
+                        yt[0] = 1; yt[1] = 0;
+                        xt[0] = -eq[1]/eq[0]; xt[1] = eq[2]/eq[0];
+                    }
+
+                    /* Find a linear representation of r using eqs[0]. */
+                    rt[0] = -(eqs[0][0]*xt[0] + eqs[0][1]*yt[0])/eqs[0][2];
+                    rt[1] = (eqs[0][3] - eqs[0][0]*xt[1] -
+                             eqs[0][1]*yt[1])/eqs[0][2];
+
+                    /* Construct the quadratic equation. */
+                    q[0] = -rt[0]*rt[0];
+                    q[1] = -2*rt[0]*rt[1];
+                    q[2] = -rt[1]*rt[1];
+                    q[0] += xt[0]*xt[0];
+                    q[1] += 2*xt[0]*(xt[1]-dots[0]->x);
+                    q[2] += (xt[1]-dots[0]->x)*(xt[1]-dots[0]->x);
+                    q[0] += yt[0]*yt[0];
+                    q[1] += 2*yt[0]*(yt[1]-dots[0]->y);
+                    q[2] += (yt[1]-dots[0]->y)*(yt[1]-dots[0]->y);
+
+                    /* And solve it. */
+                    disc = q[1]*q[1] - 4*q[0]*q[2];
+                    if (disc >= 0) {
+                        double t;
+
+                        disc = sqrt(disc);
+
+                        t = (-q[1] + disc) / (2*q[0]);
+                        cx[cn] = xt[0]*t + xt[1];
+                        cy[cn] = yt[0]*t + yt[1];
+                        cn++;
+
+                        t = (-q[1] - disc) / (2*q[0]);
+                        cx[cn] = xt[0]*t + xt[1];
+                        cy[cn] = yt[0]*t + yt[1];
+                        cn++;
+                    }
+                } else if (nedges == 1) {
+                    /*
+                     * Two dots and an edge. This one's another
+                     * quadratic equation.
+                     *
+                     * The point we want must lie on the perpendicular
+                     * bisector of the two dots; that much is obvious.
+                     * So we can construct a parametrisation of that
+                     * bisecting line, giving linear formulae for x,y
+                     * in terms of t. We can also express the distance
+                     * from the edge as such a linear formula.
+                     *
+                     * Then we set that equal to the radius of the
+                     * circle passing through the two points, which is
+                     * a Pythagoras exercise; that gives rise to a
+                     * quadratic in t, which we solve.
+                     */
+                    double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
+                    double q[3];                /* a,b,c: at^2+bt+c=0 */
+                    double disc;
+                    double halfsep;
+
+                    /* Find parametric formulae for x,y. */
+                    {
+                        int x1 = dots[0]->x, x2 = dots[1]->x;
+                        int y1 = dots[0]->y, y2 = dots[1]->y;
+                        int dx = x2-x1, dy = y2-y1;
+                        double d = sqrt((double)dx*dx + (double)dy*dy);
+
+                        xt[1] = (x1+x2)/2.0;
+                        yt[1] = (y1+y2)/2.0;
+                        /* It's convenient if we have t at standard scale. */
+                        xt[0] = -dy/d;
+                        yt[0] = dx/d;
+
+                        /* Also note down half the separation between
+                         * the dots, for use in computing the circle radius. */
+                        halfsep = 0.5*d;
+                    }
+
+                    /* Find a parametric formula for r. */
+                    {
+                        int x1 = edgedot1[0]->x, x2 = edgedot2[0]->x;
+                        int y1 = edgedot1[0]->y, y2 = edgedot2[0]->y;
+                        int dx = x2-x1, dy = y2-y1;
+                        double d = sqrt((double)dx*dx + (double)dy*dy);
+                        rt[0] = (xt[0]*dy - yt[0]*dx) / d;
+                        rt[1] = ((xt[1]-x1)*dy - (yt[1]-y1)*dx) / d;
+                    }
+
+                    /* Construct the quadratic equation. */
+                    q[0] = rt[0]*rt[0];
+                    q[1] = 2*rt[0]*rt[1];
+                    q[2] = rt[1]*rt[1];
+                    q[0] -= 1;
+                    q[2] -= halfsep*halfsep;
+
+                    /* And solve it. */
+                    disc = q[1]*q[1] - 4*q[0]*q[2];
+                    if (disc >= 0) {
+                        double t;
+
+                        disc = sqrt(disc);
+
+                        t = (-q[1] + disc) / (2*q[0]);
+                        cx[cn] = xt[0]*t + xt[1];
+                        cy[cn] = yt[0]*t + yt[1];
+                        cn++;
+
+                        t = (-q[1] - disc) / (2*q[0]);
+                        cx[cn] = xt[0]*t + xt[1];
+                        cy[cn] = yt[0]*t + yt[1];
+                        cn++;
+                    }
+                } else if (nedges == 0) {
+                    /*
+                     * Three dots. This is another linear matrix
+                     * equation, this time with each row of the matrix
+                     * representing the perpendicular bisector between
+                     * two of the points. Of course we only need two
+                     * such lines to find their intersection, so we
+                     * need only solve a 2x2 matrix equation.
+                     */
+
+                    double matrix[4], vector[2], vector2[2];
+                    int m;
+
+                    for (m = 0; m < 2; m++) {
+                        int x1 = dots[m]->x, x2 = dots[m+1]->x;
+                        int y1 = dots[m]->y, y2 = dots[m+1]->y;
+                        int dx = x2-x1, dy = y2-y1;
+
+                        /*
+                         * ((x,y) - (x1,y1)) . (dx,dy) = 1/2 |(dx,dy)|^2
+                         *
+                         * => 2x dx + 2y dy = dx^2+dy^2 + (2 x1 dx + 2 y1 dy)
+                         */
+                        matrix[2*m+0] = 2*dx;
+                        matrix[2*m+1] = 2*dy;
+                        vector[m] = ((double)dx*dx + (double)dy*dy +
+                                     2.0*x1*dx + 2.0*y1*dy);
+                    }
+
+                    if (solve_2x2_matrix(matrix, vector, vector2)) {
+                        cx[cn] = vector2[0];
+                        cy[cn] = vector2[1];
+                        cn++;
+                    }
+                }
+
+                /*
+                 * Now go through our candidate points and see if any
+                 * of them are better than what we've got so far.
+                 */
+                for (m = 0; m < cn; m++) {
+                    double x = cx[m], y = cy[m];
+
+                    /*
+                     * First, disqualify the point if it's not inside
+                     * the polygon, which we work out by counting the
+                     * edges to the right of the point. (For
+                     * tiebreaking purposes when edges start or end on
+                     * our y-coordinate or go right through it, we
+                     * consider our point to be offset by a small
+                     * _positive_ epsilon in both the x- and
+                     * y-direction.)
+                     */
+                    int e, in = 0;
+                    for (e = 0; e < f->order; e++) {
+                        int xs = f->edges[e]->dot1->x;
+                        int xe = f->edges[e]->dot2->x;
+                        int ys = f->edges[e]->dot1->y;
+                        int ye = f->edges[e]->dot2->y;
+                        if ((y >= ys && y < ye) || (y >= ye && y < ys)) {
+                            /*
+                             * The line goes past our y-position. Now we need
+                             * to know if its x-coordinate when it does so is
+                             * to our right.
+                             *
+                             * The x-coordinate in question is mathematically
+                             * (y - ys) * (xe - xs) / (ye - ys), and we want
+                             * to know whether (x - xs) >= that. Of course we
+                             * avoid the division, so we can work in integers;
+                             * to do this we must multiply both sides of the
+                             * inequality by ye - ys, which means we must
+                             * first check that's not negative.
+                             */
+                            int num = xe - xs, denom = ye - ys;
+                            if (denom < 0) {
+                                num = -num;
+                                denom = -denom;
+                            }
+                            if ((x - xs) * denom >= (y - ys) * num)
+                                in ^= 1;
+                        }
+                    }
+
+                    if (in) {
+#ifdef HUGE_VAL
+                        double mindist = HUGE_VAL;
+#else
+#ifdef DBL_MAX
+                        double mindist = DBL_MAX;
+#else
+#error No way to get maximum floating-point number.
+#endif
+#endif
+                        int e, d;
+
+                        /*
+                         * This point is inside the polygon, so now we check
+                         * its minimum distance to every edge and corner.
+                         * First the corners ...
+                         */
+                        for (d = 0; d < f->order; d++) {
+                            int xp = f->dots[d]->x;
+                            int yp = f->dots[d]->y;
+                            double dx = x - xp, dy = y - yp;
+                            double dist = dx*dx + dy*dy;
+                            if (mindist > dist)
+                                mindist = dist;
+                        }
+
+                        /*
+                         * ... and now also check the perpendicular distance
+                         * to every edge, if the perpendicular lies between
+                         * the edge's endpoints.
+                         */
+                        for (e = 0; e < f->order; e++) {
+                            int xs = f->edges[e]->dot1->x;
+                            int xe = f->edges[e]->dot2->x;
+                            int ys = f->edges[e]->dot1->y;
+                            int ye = f->edges[e]->dot2->y;
+
+                            /*
+                             * If s and e are our endpoints, and p our
+                             * candidate circle centre, the foot of a
+                             * perpendicular from p to the line se lies
+                             * between s and e if and only if (p-s).(e-s) lies
+                             * strictly between 0 and (e-s).(e-s).
+                             */
+                            int edx = xe - xs, edy = ye - ys;
+                            double pdx = x - xs, pdy = y - ys;
+                            double pde = pdx * edx + pdy * edy;
+                            long ede = (long)edx * edx + (long)edy * edy;
+                            if (0 < pde && pde < ede) {
+                                /*
+                                 * Yes, the nearest point on this edge is
+                                 * closer than either endpoint, so we must
+                                 * take it into account by measuring the
+                                 * perpendicular distance to the edge and
+                                 * checking its square against mindist.
+                                 */
+
+                                double pdre = pdx * edy - pdy * edx;
+                                double sqlen = pdre * pdre / ede;
+
+                                if (mindist > sqlen)
+                                    mindist = sqlen;
+                            }
+                        }
+
+                        /*
+                         * Right. Now we know the biggest circle around this
+                         * point, so we can check it against bestdist.
+                         */
+                        if (bestdist < mindist) {
+                            bestdist = mindist;
+                            xbest = x;
+                            ybest = y;
+                        }
+                    }
+                }
+
+                if (k < f->order)
+                    nedges--;
+                else
+                    ndots--;
+            }
+            if (j < f->order)
+                nedges--;
+            else
+                ndots--;
+        }
+        if (i < f->order)
+            nedges--;
+        else
+            ndots--;
+    }
+
+    assert(bestdist > 0);
+
+    f->has_incentre = TRUE;
+    f->ix = xbest + 0.5;               /* round to nearest */
+    f->iy = ybest + 0.5;
+}
+
 /* ------ Generate various types of grid ------ */
 
 /* General method is to generate faces, by calculating their dot coordinates.
@@ -682,11 +1385,23 @@ static void grid_face_set_dot(grid *g, grid_dot *d, int position)
  * arithmetic here!
  */
 
-grid *grid_new_square(int width, int height)
+#define SQUARE_TILESIZE 20
+
+static void grid_size_square(int width, int height,
+                      int *tilesize, int *xextent, int *yextent)
+{
+    int a = SQUARE_TILESIZE;
+
+    *tilesize = a;
+    *xextent = width * a;
+    *yextent = height * a;
+}
+
+static grid *grid_new_square(int width, int height, const char *desc)
 {
     int x, y;
     /* Side length */
-    int a = 20;
+    int a = SQUARE_TILESIZE;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = width * height;
@@ -694,7 +1409,7 @@ grid *grid_new_square(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
+    grid *g = grid_empty();
     g->tilesize = a;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
@@ -724,27 +1439,41 @@ grid *grid_new_square(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_honeycomb(int width, int height)
+#define HONEY_TILESIZE 45
+/* Vector for side of hexagon - ratio is close to sqrt(3) */
+#define HONEY_A 15
+#define HONEY_B 26
+
+static void grid_size_honeycomb(int width, int height,
+                         int *tilesize, int *xextent, int *yextent)
+{
+    int a = HONEY_A;
+    int b = HONEY_B;
+
+    *tilesize = HONEY_TILESIZE;
+    *xextent = (3 * a * (width-1)) + 4*a;
+    *yextent = (2 * b * (height-1)) + 3*b;
+}
+
+static grid *grid_new_honeycomb(int width, int height, const char *desc)
 {
     int x, y;
-    /* Vector for side of hexagon - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = HONEY_A;
+    int b = HONEY_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = width * height;
     int max_dots = 2 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 3 * a;
+    grid *g = grid_empty();
+    g->tilesize = HONEY_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -779,108 +1508,238 @@ grid *grid_new_honeycomb(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
+#define TRIANGLE_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define TRIANGLE_VEC_X 15
+#define TRIANGLE_VEC_Y 26
+
+static void grid_size_triangular(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int vec_x = TRIANGLE_VEC_X;
+    int vec_y = TRIANGLE_VEC_Y;
+
+    *tilesize = TRIANGLE_TILESIZE;
+    *xextent = (width+1) * 2 * vec_x;
+    *yextent = height * vec_y;
+}
+
+static char *grid_validate_desc_triangular(grid_type type, int width,
+                                           int height, const char *desc)
+{
+    /*
+     * Triangular grids: an absent description is valid (indicating
+     * the old-style approach which had 'ears', i.e. triangles
+     * connected to only one other face, at some grid corners), and so
+     * is a description reading just "0" (indicating the new-style
+     * approach in which those ears are trimmed off). Anything else is
+     * illegal.
+     */
+
+    if (!desc || !strcmp(desc, "0"))
+        return NULL;
+
+    return "Unrecognised grid description.";
+}
+
 /* Doesn't use the previous method of generation, it pre-dates it!
  * A triangular grid is just about simple enough to do by "brute force" */
-grid *grid_new_triangular(int width, int height)
+static grid *grid_new_triangular(int width, int height, const char *desc)
 {
     int x,y;
+    int version = (desc == NULL ? -1 : atoi(desc));
     
     /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int vec_x = 15;
-    int vec_y = 26;
+    int vec_x = TRIANGLE_VEC_X;
+    int vec_y = TRIANGLE_VEC_Y;
     
     int index;
 
     /* convenient alias */
     int w = width + 1;
 
-    grid *g = grid_new();
-    g->tilesize = 18; /* adjust to your taste */
-
-    g->num_faces = width * height * 2;
-    g->num_dots = (width + 1) * (height + 1);
-    g->faces = snewn(g->num_faces, grid_face);
-    g->dots = snewn(g->num_dots, grid_dot);
-
-    /* generate dots */
-    index = 0;
-    for (y = 0; y <= height; y++) {
-        for (x = 0; x <= width; x++) {
-            grid_dot *d = g->dots + index;
-            /* odd rows are offset to the right */
-            d->order = 0;
-            d->edges = NULL;
-            d->faces = NULL;
-            d->x = x * 2 * vec_x + ((y % 2) ? vec_x : 0);
-            d->y = y * vec_y;
-            index++;
+    grid *g = grid_empty();
+    g->tilesize = TRIANGLE_TILESIZE;
+
+    if (version == -1) {
+        /*
+         * Old-style triangular grid generation, preserved as-is for
+         * backwards compatibility with old game ids, in which it's
+         * just a little asymmetric and there are 'ears' (faces linked
+         * to only one other face) at two grid corners.
+         *
+         * Example old-style game ids, which should still work, and in
+         * which you should see the ears in the TL/BR corners on the
+         * first one and in the TL/BL corners on the second:
+         *
+         *   5x5t1:2c2a1a2a201a1a1a1112a1a2b1211f0b21a2a2a0a
+         *   5x6t1:a022a212h1a1d1a12c2b11a012b1a20d1a0a12e
+         */
+
+        g->num_faces = width * height * 2;
+        g->num_dots = (width + 1) * (height + 1);
+        g->faces = snewn(g->num_faces, grid_face);
+        g->dots = snewn(g->num_dots, grid_dot);
+
+        /* generate dots */
+        index = 0;
+        for (y = 0; y <= height; y++) {
+            for (x = 0; x <= width; x++) {
+                grid_dot *d = g->dots + index;
+                /* odd rows are offset to the right */
+                d->order = 0;
+                d->edges = NULL;
+                d->faces = NULL;
+                d->x = x * 2 * vec_x + ((y % 2) ? vec_x : 0);
+                d->y = y * vec_y;
+                index++;
+            }
         }
-    }
     
-    /* generate faces */
-    index = 0;
-    for (y = 0; y < height; y++) {
-        for (x = 0; x < width; x++) {
-            /* initialise two faces for this (x,y) */
-            grid_face *f1 = g->faces + index;
-            grid_face *f2 = f1 + 1;
-            f1->edges = NULL;
-            f1->order = 3;
-            f1->dots = snewn(f1->order, grid_dot*);
-            f2->edges = NULL;
-            f2->order = 3;
-            f2->dots = snewn(f2->order, grid_dot*);
-
-            /* face descriptions depend on whether the row-number is
-             * odd or even */
+        /* generate faces */
+        index = 0;
+        for (y = 0; y < height; y++) {
+            for (x = 0; x < width; x++) {
+                /* initialise two faces for this (x,y) */
+                grid_face *f1 = g->faces + index;
+                grid_face *f2 = f1 + 1;
+                f1->edges = NULL;
+                f1->order = 3;
+                f1->dots = snewn(f1->order, grid_dot*);
+                f1->has_incentre = FALSE;
+                f2->edges = NULL;
+                f2->order = 3;
+                f2->dots = snewn(f2->order, grid_dot*);
+                f2->has_incentre = FALSE;
+
+                /* face descriptions depend on whether the row-number is
+                 * odd or even */
+                if (y % 2) {
+                    f1->dots[0] = g->dots + y       * w + x;
+                    f1->dots[1] = g->dots + (y + 1) * w + x + 1;
+                    f1->dots[2] = g->dots + (y + 1) * w + x;
+                    f2->dots[0] = g->dots + y       * w + x;
+                    f2->dots[1] = g->dots + y       * w + x + 1;
+                    f2->dots[2] = g->dots + (y + 1) * w + x + 1;
+                } else {
+                    f1->dots[0] = g->dots + y       * w + x;
+                    f1->dots[1] = g->dots + y       * w + x + 1;
+                    f1->dots[2] = g->dots + (y + 1) * w + x;
+                    f2->dots[0] = g->dots + y       * w + x + 1;
+                    f2->dots[1] = g->dots + (y + 1) * w + x + 1;
+                    f2->dots[2] = g->dots + (y + 1) * w + x;
+                }
+                index += 2;
+            }
+        }
+    } else {
+        /*
+         * New-style approach, in which there are never any 'ears',
+         * and if height is even then the grid is nicely 4-way
+         * symmetric.
+         *
+         * Example new-style grids:
+         *
+         *   5x5t1:0_21120b11a1a01a1a00c1a0b211021c1h1a2a1a0a
+         *   5x6t1:0_a1212c22c2a02a2f22a0c12a110d0e1c0c0a101121a1
+         */
+        tree234 *points = newtree234(grid_point_cmp_fn);
+        /* Upper bounds - don't have to be exact */
+        int max_faces = height * (2*width+1);
+        int max_dots = (height+1) * (width+1) * 4;
+
+        g->faces = snewn(max_faces, grid_face);
+        g->dots = snewn(max_dots, grid_dot);
+
+        for (y = 0; y < height; y++) {
+            /*
+             * Each row contains (width+1) triangles one way up, and
+             * (width) triangles the other way up. Which way up is
+             * which varies with parity of y. Also, the dots around
+             * each face will flip direction with parity of y, so we
+             * set up n1 and n2 to cope with that easily.
+             */
+            int y0, y1, n1, n2;
+            y0 = y1 = y * vec_y;
             if (y % 2) {
-                f1->dots[0] = g->dots + y       * w + x;
-                f1->dots[1] = g->dots + (y + 1) * w + x + 1;
-                f1->dots[2] = g->dots + (y + 1) * w + x;
-                f2->dots[0] = g->dots + y       * w + x;
-                f2->dots[1] = g->dots + y       * w + x + 1;
-                f2->dots[2] = g->dots + (y + 1) * w + x + 1;
+                y1 += vec_y;
+                n1 = 2; n2 = 1;
             } else {
-                f1->dots[0] = g->dots + y       * w + x;
-                f1->dots[1] = g->dots + y       * w + x + 1;
-                f1->dots[2] = g->dots + (y + 1) * w + x;
-                f2->dots[0] = g->dots + y       * w + x + 1;
-                f2->dots[1] = g->dots + (y + 1) * w + x + 1;
-                f2->dots[2] = g->dots + (y + 1) * w + x;
+                y0 += vec_y;
+                n1 = 1; n2 = 2;
+            }
+
+            for (x = 0; x <= width; x++) {
+                int x0 = 2*x * vec_x, x1 = x0 + vec_x, x2 = x1 + vec_x;
+
+                /*
+                 * If the grid has odd height, then we skip the first
+                 * and last triangles on this row, otherwise they'll
+                 * end up as ears.
+                 */
+                if (height % 2 == 1 && y == height-1 && (x == 0 || x == width))
+                    continue;
+
+                grid_face_add_new(g, 3);
+                grid_face_set_dot(g, grid_get_dot(g, points, x0, y0), 0);
+                grid_face_set_dot(g, grid_get_dot(g, points, x1, y1), n1);
+                grid_face_set_dot(g, grid_get_dot(g, points, x2, y0), n2);
+            }
+
+            for (x = 0; x < width; x++) {
+                int x0 = (2*x+1) * vec_x, x1 = x0 + vec_x, x2 = x1 + vec_x;
+
+                grid_face_add_new(g, 3);
+                grid_face_set_dot(g, grid_get_dot(g, points, x0, y1), 0);
+                grid_face_set_dot(g, grid_get_dot(g, points, x1, y0), n2);
+                grid_face_set_dot(g, grid_get_dot(g, points, x2, y1), n1);
             }
-            index += 2;
         }
-    }
 
-    /* "+ width" takes us to the middle of the row, because each row has
-     * (2*width) faces. */
-    g->middle_face = g->faces + (height / 2) * 2 * width + width;
+        freetree234(points);
+        assert(g->num_faces <= max_faces);
+        assert(g->num_dots <= max_dots);
+    }
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_snubsquare(int width, int height)
+#define SNUBSQUARE_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define SNUBSQUARE_A 15
+#define SNUBSQUARE_B 26
+
+static void grid_size_snubsquare(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = SNUBSQUARE_A;
+    int b = SNUBSQUARE_B;
+
+    *tilesize = SNUBSQUARE_TILESIZE;
+    *xextent = (a+b) * (width-1) + a + b;
+    *yextent = (a+b) * (height-1) + a + b;
+}
+
+static grid *grid_new_snubsquare(int width, int height, const char *desc)
 {
     int x, y;
-    /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = SNUBSQUARE_A;
+    int b = SNUBSQUARE_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 3 * width * height;
     int max_dots = 2 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 18;
+    grid *g = grid_empty();
+    g->tilesize = SNUBSQUARE_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -960,27 +1819,40 @@ grid *grid_new_snubsquare(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_cairo(int width, int height)
+#define CAIRO_TILESIZE 40
+/* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
+#define CAIRO_A 14
+#define CAIRO_B 31
+
+static void grid_size_cairo(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int b = CAIRO_B; /* a unused in determining grid size. */
+
+    *tilesize = CAIRO_TILESIZE;
+    *xextent = 2*b*(width-1) + 2*b;
+    *yextent = 2*b*(height-1) + 2*b;
+}
+
+static grid *grid_new_cairo(int width, int height, const char *desc)
 {
     int x, y;
-    /* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
-    int a = 14;
-    int b = 31;
+    int a = CAIRO_A;
+    int b = CAIRO_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 2 * width * height;
     int max_dots = 3 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = CAIRO_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1053,18 +1925,32 @@ grid *grid_new_cairo(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_greathexagonal(int width, int height)
+#define GREATHEX_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define GREATHEX_A 15
+#define GREATHEX_B 26
+
+static void grid_size_greathexagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = GREATHEX_A;
+    int b = GREATHEX_B;
+
+    *tilesize = GREATHEX_TILESIZE;
+    *xextent = (3*a + b) * (width-1) + 4*a;
+    *yextent = (2*a + 2*b) * (height-1) + 3*b + a;
+}
+
+static grid *grid_new_greathexagonal(int width, int height, const char *desc)
 {
     int x, y;
-    /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = GREATHEX_A;
+    int b = GREATHEX_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 6 * (width + 1) * (height + 1);
@@ -1072,8 +1958,8 @@ grid *grid_new_greathexagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 18;
+    grid *g = grid_empty();
+    g->tilesize = GREATHEX_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1169,18 +2055,32 @@ grid *grid_new_greathexagonal(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_octagonal(int width, int height)
+#define OCTAGONAL_TILESIZE 40
+/* b/a approx sqrt(2) */
+#define OCTAGONAL_A 29
+#define OCTAGONAL_B 41
+
+static void grid_size_octagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = OCTAGONAL_A;
+    int b = OCTAGONAL_B;
+
+    *tilesize = OCTAGONAL_TILESIZE;
+    *xextent = (2*a + b) * width;
+    *yextent = (2*a + b) * height;
+}
+
+static grid *grid_new_octagonal(int width, int height, const char *desc)
 {
     int x, y;
-    /* b/a approx sqrt(2) */
-    int a = 29;
-    int b = 41;
+    int a = OCTAGONAL_A;
+    int b = OCTAGONAL_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 2 * width * height;
@@ -1188,8 +2088,8 @@ grid *grid_new_octagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = OCTAGONAL_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1238,18 +2138,32 @@ grid *grid_new_octagonal(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + (height/2) * width + (width/2);
 
     grid_make_consistent(g);
     return g;
 }
 
-grid *grid_new_kites(int width, int height)
+#define KITE_TILESIZE 40
+/* b/a approx sqrt(3) */
+#define KITE_A 15
+#define KITE_B 26
+
+static void grid_size_kites(int width, int height,
+                     int *tilesize, int *xextent, int *yextent)
+{
+    int a = KITE_A;
+    int b = KITE_B;
+
+    *tilesize = KITE_TILESIZE;
+    *xextent = 4*b * width + 2*b;
+    *yextent = 6*a * (height-1) + 8*a;
+}
+
+static grid *grid_new_kites(int width, int height, const char *desc)
 {
     int x, y;
-    /* b/a approx sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = KITE_A;
+    int b = KITE_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 6 * width * height;
@@ -1257,8 +2171,8 @@ grid *grid_new_kites(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = KITE_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1344,10 +2258,639 @@ grid *grid_new_kites(int width, int height)
     freetree234(points);
     assert(g->num_faces <= max_faces);
     assert(g->num_dots <= max_dots);
-    g->middle_face = g->faces + 6 * ((height/2) * width + (width/2));
 
     grid_make_consistent(g);
     return g;
 }
 
+#define FLORET_TILESIZE 150
+/* -py/px is close to tan(30 - atan(sqrt(3)/9))
+ * using py=26 makes everything lean to the left, rather than right
+ */
+#define FLORET_PX 75
+#define FLORET_PY -26
+
+static void grid_size_floret(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
+    int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
+    int ry = qy-py;
+    /* rx unused in determining grid size. */
+
+    *tilesize = FLORET_TILESIZE;
+    *xextent = (6*px+3*qx)/2 * (width-1) + 4*qx + 2*px;
+    *yextent = (5*qy-4*py) * (height-1) + 4*qy + 2*ry;
+}
+
+static grid *grid_new_floret(int width, int height, const char *desc)
+{
+    int x, y;
+    /* Vectors for sides; weird numbers needed to keep puzzle aligned with window
+     * -py/px is close to tan(30 - atan(sqrt(3)/9))
+     * using py=26 makes everything lean to the left, rather than right
+     */
+    int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
+    int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
+    int rx = qx-px, ry = qy-py;                 /* |(-15,  78)| = 79.38 */
+
+    /* Upper bounds - don't have to be exact */
+    int max_faces = 6 * width * height;
+    int max_dots = 9 * (width + 1) * (height + 1);
+
+    tree234 *points;
+
+    grid *g = grid_empty();
+    g->tilesize = FLORET_TILESIZE;
+    g->faces = snewn(max_faces, grid_face);
+    g->dots = snewn(max_dots, grid_dot);
+
+    points = newtree234(grid_point_cmp_fn);
+
+    /* generate pentagonal faces */
+    for (y = 0; y < height; y++) {
+        for (x = 0; x < width; x++) {
+            grid_dot *d;
+            /* face centre */
+            int cx = (6*px+3*qx)/2 * x;
+            int cy = (4*py-5*qy) * y;
+            if (x % 2)
+                cy -= (4*py-5*qy)/2;
+            else if (y && y == height-1)
+                continue; /* make better looking grids?  try 3x3 for instance */
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx+2*rx+qx, cy+2*ry+qy); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx+2*qx+rx, cy+2*qy+ry); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 4);
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx+2*qx+px, cy+2*qy+py); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx+2*px+qx, cy+2*py+qy); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 4);
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx+2*px-rx, cy+2*py-ry); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx-2*rx+px, cy-2*ry+py); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 4);
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx-2*rx-qx, cy-2*ry-qy); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx-2*qx-rx, cy-2*qy-ry); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 4);
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx-2*qx-px, cy-2*qy-py); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx-2*px-qx, cy-2*py-qy); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 4);
+
+            grid_face_add_new(g, 5);
+            d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, cx-2*px+rx, cy-2*py+ry); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, cx+2*rx-px, cy+2*ry-py); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 4);
+        }
+    }
+
+    freetree234(points);
+    assert(g->num_faces <= max_faces);
+    assert(g->num_dots <= max_dots);
+
+    grid_make_consistent(g);
+    return g;
+}
+
+/* DODEC_* are used for dodecagonal and great-dodecagonal grids. */
+#define DODEC_TILESIZE 26
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define DODEC_A 15
+#define DODEC_B 26
+
+static void grid_size_dodecagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    *tilesize = DODEC_TILESIZE;
+    *xextent = (4*a + 2*b) * (width-1) + 3*(2*a + b);
+    *yextent = (3*a + 2*b) * (height-1) + 2*(2*a + b);
+}
+
+static grid *grid_new_dodecagonal(int width, int height, const char *desc)
+{
+    int x, y;
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    /* Upper bounds - don't have to be exact */
+    int max_faces = 3 * width * height;
+    int max_dots = 14 * width * height;
+
+    tree234 *points;
+
+    grid *g = grid_empty();
+    g->tilesize = DODEC_TILESIZE;
+    g->faces = snewn(max_faces, grid_face);
+    g->dots = snewn(max_dots, grid_dot);
+
+    points = newtree234(grid_point_cmp_fn);
+
+    for (y = 0; y < height; y++) {
+        for (x = 0; x < width; x++) {
+            grid_dot *d;
+            /* centre of dodecagon */
+            int px = (4*a + 2*b) * x;
+            int py = (3*a + 2*b) * y;
+            if (y % 2)
+                px += 2*a + b;
+
+            /* dodecagon */
+            grid_face_add_new(g, 12);
+            d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
+            d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
+            d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
+            d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
+            d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
+            d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
+            d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
+            d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
+
+            /* triangle below dodecagon */
+           if ((y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
+               grid_face_add_new(g, 3);
+               d = grid_get_dot(g, points, px + a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px    , py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px - a, py + (2*a +   b)); grid_face_set_dot(g, d, 2);
+           }
+
+            /* triangle above dodecagon */
+           if ((y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
+               grid_face_add_new(g, 3);
+               d = grid_get_dot(g, points, px - a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px    , py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px + a, py - (2*a +   b)); grid_face_set_dot(g, d, 2);
+           }
+       }
+    }
+
+    freetree234(points);
+    assert(g->num_faces <= max_faces);
+    assert(g->num_dots <= max_dots);
+
+    grid_make_consistent(g);
+    return g;
+}
+
+static void grid_size_greatdodecagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    *tilesize = DODEC_TILESIZE;
+    *xextent = (6*a + 2*b) * (width-1) + 2*(2*a + b) + 3*a + b;
+    *yextent = (3*a + 3*b) * (height-1) + 2*(2*a + b);
+}
+
+static grid *grid_new_greatdodecagonal(int width, int height, const char *desc)
+{
+    int x, y;
+    /* Vector for side of triangle - ratio is close to sqrt(3) */
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    /* Upper bounds - don't have to be exact */
+    int max_faces = 30 * width * height;
+    int max_dots = 200 * width * height;
+
+    tree234 *points;
+
+    grid *g = grid_empty();
+    g->tilesize = DODEC_TILESIZE;
+    g->faces = snewn(max_faces, grid_face);
+    g->dots = snewn(max_dots, grid_dot);
+
+    points = newtree234(grid_point_cmp_fn);
+
+    for (y = 0; y < height; y++) {
+        for (x = 0; x < width; x++) {
+            grid_dot *d;
+            /* centre of dodecagon */
+            int px = (6*a + 2*b) * x;
+            int py = (3*a + 3*b) * y;
+            if (y % 2)
+                px += 3*a + b;
+
+            /* dodecagon */
+            grid_face_add_new(g, 12);
+            d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
+            d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
+            d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
+            d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
+            d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
+            d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
+            d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
+            d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
+            d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
+            d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
+            d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
+            d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
+
+            /* hexagon below dodecagon */
+           if (y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
+               grid_face_add_new(g, 6);
+               d = grid_get_dot(g, points, px +   a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px + 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px +   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 2);
+               d = grid_get_dot(g, points, px -   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 3);
+               d = grid_get_dot(g, points, px - 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 4);
+               d = grid_get_dot(g, points, px -   a, py + (2*a +   b)); grid_face_set_dot(g, d, 5);
+           }
+
+            /* hexagon above dodecagon */
+           if (y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
+               grid_face_add_new(g, 6);
+               d = grid_get_dot(g, points, px -   a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px - 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px -   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 2);
+               d = grid_get_dot(g, points, px +   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 3);
+               d = grid_get_dot(g, points, px + 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 4);
+               d = grid_get_dot(g, points, px +   a, py - (2*a +   b)); grid_face_set_dot(g, d, 5);
+           }
+
+            /* square on right of dodecagon */
+           if (x < width - 1) {
+               grid_face_add_new(g, 4);
+               d = grid_get_dot(g, points, px + 2*a + b, py - a); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px + 4*a + b, py - a); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px + 4*a + b, py + a); grid_face_set_dot(g, d, 2);
+               d = grid_get_dot(g, points, px + 2*a + b, py + a); grid_face_set_dot(g, d, 3);
+           }
+
+            /* square on top right of dodecagon */
+           if (y && (x < width - 1 || !(y % 2))) {
+               grid_face_add_new(g, 4);
+               d = grid_get_dot(g, points, px + (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px + (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px + (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 2);
+               d = grid_get_dot(g, points, px + (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 3);
+           }
+
+            /* square on top left of dodecagon */
+           if (y && (x || (y % 2))) {
+               grid_face_add_new(g, 4);
+               d = grid_get_dot(g, points, px - (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 0);
+               d = grid_get_dot(g, points, px - (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 1);
+               d = grid_get_dot(g, points, px - (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 2);
+               d = grid_get_dot(g, points, px - (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 3);
+           }
+       }
+    }
+
+    freetree234(points);
+    assert(g->num_faces <= max_faces);
+    assert(g->num_dots <= max_dots);
+
+    grid_make_consistent(g);
+    return g;
+}
+
+typedef struct setface_ctx
+{
+    int xmin, xmax, ymin, ymax;
+
+    grid *g;
+    tree234 *points;
+} setface_ctx;
+
+static double round_int_nearest_away(double r)
+{
+    return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
+}
+
+static int set_faces(penrose_state *state, vector *vs, int n, int depth)
+{
+    setface_ctx *sf_ctx = (setface_ctx *)state->ctx;
+    int i;
+    int xs[4], ys[4];
+
+    if (depth < state->max_depth) return 0;
+#ifdef DEBUG_PENROSE
+    if (n != 4) return 0; /* triangles are sent as debugging. */
+#endif
+
+    for (i = 0; i < n; i++) {
+        double tx = v_x(vs, i), ty = v_y(vs, i);
+
+        xs[i] = (int)round_int_nearest_away(tx);
+        ys[i] = (int)round_int_nearest_away(ty);
+
+        if (xs[i] < sf_ctx->xmin || xs[i] > sf_ctx->xmax) return 0;
+        if (ys[i] < sf_ctx->ymin || ys[i] > sf_ctx->ymax) return 0;
+    }
+
+    grid_face_add_new(sf_ctx->g, n);
+    debug(("penrose: new face l=%f gen=%d...",
+           penrose_side_length(state->start_size, depth), depth));
+    for (i = 0; i < n; i++) {
+        grid_dot *d = grid_get_dot(sf_ctx->g, sf_ctx->points,
+                                   xs[i], ys[i]);
+        grid_face_set_dot(sf_ctx->g, d, i);
+        debug((" ... dot 0x%x (%d,%d) (was %2.2f,%2.2f)",
+               d, d->x, d->y, v_x(vs, i), v_y(vs, i)));
+    }
+
+    return 0;
+}
+
+#define PENROSE_TILESIZE 100
+
+static void grid_size_penrose(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    int l = PENROSE_TILESIZE;
+
+    *tilesize = l;
+    *xextent = l * width;
+    *yextent = l * height;
+}
+
+static grid *grid_new_penrose(int width, int height, int which, const char *desc); /* forward reference */
+
+static char *grid_new_desc_penrose(grid_type type, int width, int height, random_state *rs)
+{
+    int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff;
+    double outer_radius;
+    int inner_radius;
+    char gd[255];
+    int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+    grid *g;
+
+    while (1) {
+        /* We want to produce a random bit of penrose tiling, so we
+         * calculate a random offset (within the patch that penrose.c
+         * calculates for us) and an angle (multiple of 36) to rotate
+         * the patch. */
+
+        penrose_calculate_size(which, tilesize, width, height,
+                               &outer_radius, &startsz, &depth);
+
+        /* Calculate radius of (circumcircle of) patch, subtract from
+         * radius calculated. */
+        inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
+
+        /* Pick a random offset (the easy way: choose within outer
+         * square, discarding while it's outside the circle) */
+        do {
+            xoff = random_upto(rs, 2*inner_radius) - inner_radius;
+            yoff = random_upto(rs, 2*inner_radius) - inner_radius;
+        } while (sqrt(xoff*xoff+yoff*yoff) > inner_radius);
+
+        aoff = random_upto(rs, 360/36) * 36;
+
+        debug(("grid_desc: ts %d, %dx%d patch, orad %2.2f irad %d",
+               tilesize, width, height, outer_radius, inner_radius));
+        debug(("    -> xoff %d yoff %d aoff %d", xoff, yoff, aoff));
+
+        sprintf(gd, "G%d,%d,%d", xoff, yoff, aoff);
+
+        /*
+         * Now test-generate our grid, to make sure it actually
+         * produces something.
+         */
+        g = grid_new_penrose(width, height, which, gd);
+        if (g) {
+            grid_free(g);
+            break;
+        }
+        /* If not, go back to the top of this while loop and try again
+         * with a different random offset. */
+    }
+
+    return dupstr(gd);
+}
+
+static char *grid_validate_desc_penrose(grid_type type, int width, int height,
+                                        const char *desc)
+{
+    int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff, inner_radius;
+    double outer_radius;
+    int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+    grid *g;
+
+    if (!desc)
+        return "Missing grid description string.";
+
+    penrose_calculate_size(which, tilesize, width, height,
+                           &outer_radius, &startsz, &depth);
+    inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
+
+    if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
+        return "Invalid format grid description string.";
+
+    if (sqrt(xoff*xoff + yoff*yoff) > inner_radius)
+        return "Patch offset out of bounds.";
+    if ((aoff % 36) != 0 || aoff < 0 || aoff >= 360)
+        return "Angle offset out of bounds.";
+
+    /*
+     * Test-generate to ensure these parameters don't end us up with
+     * no grid at all.
+     */
+    g = grid_new_penrose(width, height, which, desc);
+    if (!g)
+        return "Patch coordinates do not identify a usable grid fragment";
+    grid_free(g);
+
+    return NULL;
+}
+
+/*
+ * We're asked for a grid of a particular size, and we generate enough
+ * of the tiling so we can be sure to have enough random grid from which
+ * to pick.
+ */
+
+static grid *grid_new_penrose(int width, int height, int which, const char *desc)
+{
+    int max_faces, max_dots, tilesize = PENROSE_TILESIZE;
+    int xsz, ysz, xoff, yoff, aoff;
+    double rradius;
+
+    tree234 *points;
+    grid *g;
+
+    penrose_state ps;
+    setface_ctx sf_ctx;
+
+    penrose_calculate_size(which, tilesize, width, height,
+                           &rradius, &ps.start_size, &ps.max_depth);
+
+    debug(("penrose: w%d h%d, tile size %d, start size %d, depth %d",
+           width, height, tilesize, ps.start_size, ps.max_depth));
+
+    ps.new_tile = set_faces;
+    ps.ctx = &sf_ctx;
+
+    max_faces = (width*3) * (height*3); /* somewhat paranoid... */
+    max_dots = max_faces * 4; /* ditto... */
+
+    g = grid_empty();
+    g->tilesize = tilesize;
+    g->faces = snewn(max_faces, grid_face);
+    g->dots = snewn(max_dots, grid_dot);
+
+    points = newtree234(grid_point_cmp_fn);
+
+    memset(&sf_ctx, 0, sizeof(sf_ctx));
+    sf_ctx.g = g;
+    sf_ctx.points = points;
+
+    if (desc != NULL) {
+        if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
+            assert(!"Invalid grid description.");
+    } else {
+        xoff = yoff = aoff = 0;
+    }
+
+    xsz = width * tilesize;
+    ysz = height * tilesize;
+
+    sf_ctx.xmin = xoff - xsz/2;
+    sf_ctx.xmax = xoff + xsz/2;
+    sf_ctx.ymin = yoff - ysz/2;
+    sf_ctx.ymax = yoff + ysz/2;
+
+    debug(("penrose: centre (%f, %f) xsz %f ysz %f",
+           0.0, 0.0, xsz, ysz));
+    debug(("penrose: x range (%f --> %f), y range (%f --> %f)",
+           sf_ctx.xmin, sf_ctx.xmax, sf_ctx.ymin, sf_ctx.ymax));
+
+    penrose(&ps, which, aoff);
+
+    freetree234(points);
+    assert(g->num_faces <= max_faces);
+    assert(g->num_dots <= max_dots);
+
+    debug(("penrose: %d faces total (equivalent to %d wide by %d high)",
+           g->num_faces, g->num_faces/height, g->num_faces/width));
+
+    /*
+     * Return NULL if we ended up with an empty grid, either because
+     * the initial generation was over too small a rectangle to
+     * encompass any face or because grid_trim_vigorously ended up
+     * removing absolutely everything.
+     */
+    if (g->num_faces == 0 || g->num_dots == 0) {
+        grid_free(g);
+        return NULL;
+    }
+    grid_trim_vigorously(g);
+    if (g->num_faces == 0 || g->num_dots == 0) {
+        grid_free(g);
+        return NULL;
+    }
+
+    grid_make_consistent(g);
+
+    /*
+     * Centre the grid in its originally promised rectangle.
+     */
+    g->lowest_x -= ((sf_ctx.xmax - sf_ctx.xmin) -
+                    (g->highest_x - g->lowest_x)) / 2;
+    g->highest_x = g->lowest_x + (sf_ctx.xmax - sf_ctx.xmin);
+    g->lowest_y -= ((sf_ctx.ymax - sf_ctx.ymin) -
+                    (g->highest_y - g->lowest_y)) / 2;
+    g->highest_y = g->lowest_y + (sf_ctx.ymax - sf_ctx.ymin);
+
+    return g;
+}
+
+static void grid_size_penrose_p2_kite(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_size_penrose(width, height, tilesize, xextent, yextent);
+}
+
+static void grid_size_penrose_p3_thick(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_size_penrose(width, height, tilesize, xextent, yextent);
+}
+
+static grid *grid_new_penrose_p2_kite(int width, int height, const char *desc)
+{
+    return grid_new_penrose(width, height, PENROSE_P2, desc);
+}
+
+static grid *grid_new_penrose_p3_thick(int width, int height, const char *desc)
+{
+    return grid_new_penrose(width, height, PENROSE_P3, desc);
+}
+
 /* ----------- End of grid generators ------------- */
+
+#define FNNEW(upper,lower) &grid_new_ ## lower,
+#define FNSZ(upper,lower) &grid_size_ ## lower,
+
+static grid *(*(grid_news[]))(int, int, const char*) = { GRIDGEN_LIST(FNNEW) };
+static void(*(grid_sizes[]))(int, int, int*, int*, int*) = { GRIDGEN_LIST(FNSZ) };
+
+char *grid_new_desc(grid_type type, int width, int height, random_state *rs)
+{
+    if (type == GRID_PENROSE_P2 || type == GRID_PENROSE_P3) {
+        return grid_new_desc_penrose(type, width, height, rs);
+    } else if (type == GRID_TRIANGULAR) {
+        return dupstr("0"); /* up-to-date version of triangular grid */
+    } else {
+        return NULL;
+    }
+}
+
+char *grid_validate_desc(grid_type type, int width, int height,
+                         const char *desc)
+{
+    if (type == GRID_PENROSE_P2 || type == GRID_PENROSE_P3) {
+        return grid_validate_desc_penrose(type, width, height, desc);
+    } else if (type == GRID_TRIANGULAR) {
+        return grid_validate_desc_triangular(type, width, height, desc);
+    } else {
+        if (desc != NULL)
+            return "Grid description strings not used with this grid type";
+        return NULL;
+    }
+}
+
+grid *grid_new(grid_type type, int width, int height, const char *desc)
+{
+    char *err = grid_validate_desc(type, width, height, desc);
+    if (err) assert(!"Invalid grid description.");
+
+    return grid_news[type](width, height, desc);
+}
+
+void grid_compute_size(grid_type type, int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_sizes[type](width, height, tilesize, xextent, yextent);
+}
+
+/* ----------- End of grid helpers ------------- */
+
+/* vim: set shiftwidth=4 tabstop=8: */