grid.c

   1 /*
   2  * (c) Lambros Lambrou 2008
   3  *
   4  * Code for working with general grids, which can be any planar graph
   5  * with faces, edges and vertices (dots).  Includes generators for a few
   6  * types of grid, including square, hexagonal, triangular and others.
   7  */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <assert.h>
  13 #include <ctype.h>
  14 #include <math.h>
  15 #include <float.h>
  16
  17 #include "puzzles.h"
  18 #include "tree234.h"
  19 #include "grid.h"
  20 #include "penrose.h"
  21
  22 /* Debugging options */
  23
  24 /*
  25 #define DEBUG_GRID
  26 */
  27
  28 /* ----------------------------------------------------------------------
  29  * Deallocate or dereference a grid
  30  */
  31 void grid_free(grid *g)
  32 {
  33     assert(g->refcount);
  34
  35     g->refcount--;
  36     if (g->refcount == 0) {
  37         int i;
  38         for (i = 0; i < g->num_faces; i++) {
  39             sfree(g->faces[i].dots);
  40             sfree(g->faces[i].edges);
  41         }
  42         for (i = 0; i < g->num_dots; i++) {
  43             sfree(g->dots[i].faces);
  44             sfree(g->dots[i].edges);
  45         }
  46         sfree(g->faces);
  47         sfree(g->edges);
  48         sfree(g->dots);
  49         sfree(g);
  50     }
  51 }
  52
  53 /* Used by the other grid generators.  Create a brand new grid with nothing
  54  * initialised (all lists are NULL) */
  55 static grid *grid_empty(void)
  56 {
  57     grid *g = snew(grid);
  58     g->faces = NULL;
  59     g->edges = NULL;
  60     g->dots = NULL;
  61     g->num_faces = g->num_edges = g->num_dots = 0;
  62     g->refcount = 1;
  63     g->lowest_x = g->lowest_y = g->highest_x = g->highest_y = 0;
  64     return g;
  65 }
  66
  67 /* Helper function to calculate perpendicular distance from
  68  * a point P to a line AB.  A and B mustn't be equal here.
  69  *
  70  * Well-known formula for area A of a triangle:
  71  *                             /  1   1   1 \
  72  * 2A = determinant of matrix  | px  ax  bx |
  73  *                             \ py  ay  by /
  74  *
  75  * Also well-known: 2A = base * height
  76  *                     = perpendicular distance * line-length.
  77  *
  78  * Combining gives: distance = determinant / line-length(a,b)
  79  */
  80 static double point_line_distance(long px, long py,
  81                                   long ax, long ay,
  82                                   long bx, long by)
  83 {
  84     long det = ax*by - bx*ay + bx*py - px*by + px*ay - ax*py;
  85     double len;
  86     det = max(det, -det);
  87     len = sqrt(SQ(ax - bx) + SQ(ay - by));
  88     return det / len;
  89 }
  90
  91 /* Determine nearest edge to where the user clicked.
  92  * (x, y) is the clicked location, converted to grid coordinates.
  93  * Returns the nearest edge, or NULL if no edge is reasonably
  94  * near the position.
  95  *
  96  * Just judging edges by perpendicular distance is not quite right -
  97  * the edge might be "off to one side". So we insist that the triangle
  98  * with (x,y) has acute angles at the edge's dots.
  99  *
 100  *     edge1
 101  *  *---------*------
 102  *            |
 103  *            |      *(x,y)
 104  *      edge2 |
 105  *            |   edge2 is OK, but edge1 is not, even though
 106  *            |   edge1 is perpendicularly closer to (x,y)
 107  *            *
 108  *
 109  */
 110 grid_edge *grid_nearest_edge(grid *g, int x, int y)
 111 {
 112     grid_edge *best_edge;
 113     double best_distance = 0;
 114     int i;
 115
 116     best_edge = NULL;
 117
 118     for (i = 0; i < g->num_edges; i++) {
 119         grid_edge *e = &g->edges[i];
 120         long e2; /* squared length of edge */
 121         long a2, b2; /* squared lengths of other sides */
 122         double dist;
 123
 124         /* See if edge e is eligible - the triangle must have acute angles
 125          * at the edge's dots.
 126          * Pythagoras formula h^2 = a^2 + b^2 detects right-angles,
 127          * so detect acute angles by testing for h^2 < a^2 + b^2 */
 128         e2 = SQ((long)e->dot1->x - (long)e->dot2->x) + SQ((long)e->dot1->y - (long)e->dot2->y);
 129         a2 = SQ((long)e->dot1->x - (long)x) + SQ((long)e->dot1->y - (long)y);
 130         b2 = SQ((long)e->dot2->x - (long)x) + SQ((long)e->dot2->y - (long)y);
 131         if (a2 >= e2 + b2) continue;
 132         if (b2 >= e2 + a2) continue;
 133
 134         /* e is eligible so far.  Now check the edge is reasonably close
 135          * to where the user clicked.  Don't want to toggle an edge if the
 136          * click was way off the grid.
 137          * There is room for experimentation here.  We could check the
 138          * perpendicular distance is within a certain fraction of the length
 139          * of the edge.  That amounts to testing a rectangular region around
 140          * the edge.
 141          * Alternatively, we could check that the angle at the point is obtuse.
 142          * That would amount to testing a circular region with the edge as
 143          * diameter. */
 144         dist = point_line_distance((long)x, (long)y,
 145                                    (long)e->dot1->x, (long)e->dot1->y,
 146                                    (long)e->dot2->x, (long)e->dot2->y);
 147         /* Is dist more than half edge length ? */
 148         if (4 * SQ(dist) > e2)
 149             continue;
 150
 151         if (best_edge == NULL || dist < best_distance) {
 152             best_edge = e;
 153             best_distance = dist;
 154         }
 155     }
 156     return best_edge;
 157 }
 158
 159 /* ----------------------------------------------------------------------
 160  * Grid generation
 161  */
 162
 163 #ifdef SVG_GRID
 164
 165 #define SVG_DOTS  1
 166 #define SVG_EDGES 2
 167 #define SVG_FACES 4
 168
 169 #define FACE_COLOUR "red"
 170 #define EDGE_COLOUR "blue"
 171 #define DOT_COLOUR "black"
 172
 173 static void grid_output_svg(FILE *fp, grid *g, int which)
 174 {
 175     int i, j;
 176
 177     fprintf(fp,"\
 178 <?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n\
 179 <!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 20010904//EN\"\n\
 180 \"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\
 181 \n\
 182 <svg xmlns=\"http://www.w3.org/2000/svg\"\n\
 183 xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n");
 184
 185     if (which & SVG_FACES) {
 186         fprintf(fp, "<g>\n");
 187         for (i = 0; i < g->num_faces; i++) {
 188             grid_face *f = g->faces + i;
 189             fprintf(fp, "<polygon points=\"");
 190             for (j = 0; j < f->order; j++) {
 191                 grid_dot *d = f->dots[j];
 192                 fprintf(fp, "%s%d,%d", (j == 0) ? "" : " ",
 193                         d->x, d->y);
 194             }
 195             fprintf(fp, "\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n",
 196                     FACE_COLOUR, FACE_COLOUR);
 197         }
 198         fprintf(fp, "</g>\n");
 199     }
 200     if (which & SVG_EDGES) {
 201         fprintf(fp, "<g>\n");
 202         for (i = 0; i < g->num_edges; i++) {
 203             grid_edge *e = g->edges + i;
 204             grid_dot *d1 = e->dot1, *d2 = e->dot2;
 205
 206             fprintf(fp, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" "
 207                         "style=\"stroke: %s\" />\n",
 208                         d1->x, d1->y, d2->x, d2->y, EDGE_COLOUR);
 209         }
 210         fprintf(fp, "</g>\n");
 211     }
 212
 213     if (which & SVG_DOTS) {
 214         fprintf(fp, "<g>\n");
 215         for (i = 0; i < g->num_dots; i++) {
 216             grid_dot *d = g->dots + i;
 217             fprintf(fp, "<ellipse cx=\"%d\" cy=\"%d\" rx=\"%d\" ry=\"%d\" fill=\"%s\" />",
 218                     d->x, d->y, g->tilesize/20, g->tilesize/20, DOT_COLOUR);
 219         }
 220         fprintf(fp, "</g>\n");
 221     }
 222
 223     fprintf(fp, "</svg>\n");
 224 }
 225 #endif
 226
 227 #ifdef SVG_GRID
 228 #include <errno.h>
 229
 230 static void grid_try_svg(grid *g, int which)
 231 {
 232     char *svg = getenv("PUZZLES_SVG_GRID");
 233     if (svg) {
 234         FILE *svgf = fopen(svg, "w");
 235         if (svgf) {
 236             grid_output_svg(svgf, g, which);
 237             fclose(svgf);
 238         } else {
 239             fprintf(stderr, "Unable to open file `%s': %s", svg, strerror(errno));
 240         }
 241     }
 242 }
 243 #endif
 244
 245 /* Show the basic grid information, before doing grid_make_consistent */
 246 static void grid_debug_basic(grid *g)
 247 {
 248     /* TODO: Maybe we should generate an SVG image of the dots and lines
 249      * of the grid here, before grid_make_consistent.
 250      * Would help with debugging grid generation. */
 251 #ifdef DEBUG_GRID
 252     int i;
 253     printf("--- Basic Grid Data ---\n");
 254     for (i = 0; i < g->num_faces; i++) {
 255         grid_face *f = g->faces + i;
 256         printf("Face %d: dots[", i);
 257         int j;
 258         for (j = 0; j < f->order; j++) {
 259             grid_dot *d = f->dots[j];
 260             printf("%s%d", j ? "," : "", (int)(d - g->dots));
 261         }
 262         printf("]\n");
 263     }
 264 #endif
 265 #ifdef SVG_GRID
 266     grid_try_svg(g, SVG_FACES);
 267 #endif
 268 }
 269
 270 /* Show the derived grid information, computed by grid_make_consistent */
 271 static void grid_debug_derived(grid *g)
 272 {
 273 #ifdef DEBUG_GRID
 274     /* edges */
 275     int i;
 276     printf("--- Derived Grid Data ---\n");
 277     for (i = 0; i < g->num_edges; i++) {
 278         grid_edge *e = g->edges + i;
 279         printf("Edge %d: dots[%d,%d] faces[%d,%d]\n",
 280             i, (int)(e->dot1 - g->dots), (int)(e->dot2 - g->dots),
 281             e->face1 ? (int)(e->face1 - g->faces) : -1,
 282             e->face2 ? (int)(e->face2 - g->faces) : -1);
 283     }
 284     /* faces */
 285     for (i = 0; i < g->num_faces; i++) {
 286         grid_face *f = g->faces + i;
 287         int j;
 288         printf("Face %d: faces[", i);
 289         for (j = 0; j < f->order; j++) {
 290             grid_edge *e = f->edges[j];
 291             grid_face *f2 = (e->face1 == f) ? e->face2 : e->face1;
 292             printf("%s%d", j ? "," : "", f2 ? (int)(f2 - g->faces) : -1);
 293         }
 294         printf("]\n");
 295     }
 296     /* dots */
 297     for (i = 0; i < g->num_dots; i++) {
 298         grid_dot *d = g->dots + i;
 299         int j;
 300         printf("Dot %d: dots[", i);
 301         for (j = 0; j < d->order; j++) {
 302             grid_edge *e = d->edges[j];
 303             grid_dot *d2 = (e->dot1 == d) ? e->dot2 : e->dot1;
 304             printf("%s%d", j ? "," : "", (int)(d2 - g->dots));
 305         }
 306         printf("] faces[");
 307         for (j = 0; j < d->order; j++) {
 308             grid_face *f = d->faces[j];
 309             printf("%s%d", j ? "," : "", f ? (int)(f - g->faces) : -1);
 310         }
 311         printf("]\n");
 312     }
 313 #endif
 314 #ifdef SVG_GRID
 315     grid_try_svg(g, SVG_DOTS | SVG_EDGES | SVG_FACES);
 316 #endif
 317 }
 318
 319 /* Helper function for building incomplete-edges list in
 320  * grid_make_consistent() */
 321 static int grid_edge_bydots_cmpfn(void *v1, void *v2)
 322 {
 323     grid_edge *a = v1;
 324     grid_edge *b = v2;
 325     grid_dot *da, *db;
 326
 327     /* Pointer subtraction is valid here, because all dots point into the
 328      * same dot-list (g->dots).
 329      * Edges are not "normalised" - the 2 dots could be stored in any order,
 330      * so we need to take this into account when comparing edges. */
 331
 332     /* Compare first dots */
 333     da = (a->dot1 < a->dot2) ? a->dot1 : a->dot2;
 334     db = (b->dot1 < b->dot2) ? b->dot1 : b->dot2;
 335     if (da != db)
 336         return db - da;
 337     /* Compare last dots */
 338     da = (a->dot1 < a->dot2) ? a->dot2 : a->dot1;
 339     db = (b->dot1 < b->dot2) ? b->dot2 : b->dot1;
 340     if (da != db)
 341         return db - da;
 342
 343     return 0;
 344 }
 345
 346 /*
 347  * 'Vigorously trim' a grid, by which I mean deleting any isolated or
 348  * uninteresting faces. By which, in turn, I mean: ensure that the
 349  * grid is composed solely of faces adjacent to at least one
 350  * 'landlocked' dot (i.e. one not in contact with the infinite
 351  * exterior face), and that all those dots are in a single connected
 352  * component.
 353  *
 354  * This function operates on, and returns, a grid satisfying the
 355  * preconditions to grid_make_consistent() rather than the
 356  * postconditions. (So call it first.)
 357  */
 358 static void grid_trim_vigorously(grid *g)
 359 {
 360     int *dotpairs, *faces, *dots;
 361     int *dsf;
 362     int i, j, k, size, newfaces, newdots;
 363
 364     /*
 365      * First construct a matrix in which each ordered pair of dots is
 366      * mapped to the index of the face in which those dots occur in
 367      * that order.
 368      */
 369     dotpairs = snewn(g->num_dots * g->num_dots, int);
 370     for (i = 0; i < g->num_dots; i++)
 371         for (j = 0; j < g->num_dots; j++)
 372             dotpairs[i*g->num_dots+j] = -1;
 373     for (i = 0; i < g->num_faces; i++) {
 374         grid_face *f = g->faces + i;
 375         int dot0 = f->dots[f->order-1] - g->dots;
 376         for (j = 0; j < f->order; j++) {
 377             int dot1 = f->dots[j] - g->dots;
 378             dotpairs[dot0 * g->num_dots + dot1] = i;
 379             dot0 = dot1;
 380         }
 381     }
 382
 383     /*
 384      * Now we can identify landlocked dots: they're the ones all of
 385      * whose edges have a mirror-image counterpart in this matrix.
 386      */
 387     dots = snewn(g->num_dots, int);
 388     for (i = 0; i < g->num_dots; i++) {
 389         dots[i] = TRUE;
 390         for (j = 0; j < g->num_dots; j++) {
 391             if ((dotpairs[i*g->num_dots+j] >= 0) ^
 392                 (dotpairs[j*g->num_dots+i] >= 0))
 393                 dots[i] = FALSE;    /* non-duplicated edge: coastal dot */
 394         }
 395     }
 396
 397     /*
 398      * Now identify connected pairs of landlocked dots, and form a dsf
 399      * unifying them.
 400      */
 401     dsf = snew_dsf(g->num_dots);
 402     for (i = 0; i < g->num_dots; i++)
 403         for (j = 0; j < i; j++)
 404             if (dots[i] && dots[j] &&
 405                 dotpairs[i*g->num_dots+j] >= 0 &&
 406                 dotpairs[j*g->num_dots+i] >= 0)
 407                 dsf_merge(dsf, i, j);
 408
 409     /*
 410      * Now look for the largest component.
 411      */
 412     size = 0;
 413     j = -1;
 414     for (i = 0; i < g->num_dots; i++) {
 415         int newsize;
 416         if (dots[i] && dsf_canonify(dsf, i) == i &&
 417             (newsize = dsf_size(dsf, i)) > size) {
 418             j = i;
 419             size = newsize;
 420         }
 421     }
 422
 423     /*
 424      * Work out which faces we're going to keep (precisely those with
 425      * at least one dot in the same connected component as j) and
 426      * which dots (those required by any face we're keeping).
 427      *
 428      * At this point we reuse the 'dots' array to indicate the dots
 429      * we're keeping, rather than the ones that are landlocked.
 430      */
 431     faces = snewn(g->num_faces, int);
 432     for (i = 0; i < g->num_faces; i++)
 433         faces[i] = 0;
 434     for (i = 0; i < g->num_dots; i++)
 435         dots[i] = 0;
 436     for (i = 0; i < g->num_faces; i++) {
 437         grid_face *f = g->faces + i;
 438         int keep = FALSE;
 439         for (k = 0; k < f->order; k++)
 440             if (dsf_canonify(dsf, f->dots[k] - g->dots) == j)
 441                 keep = TRUE;
 442         if (keep) {
 443             faces[i] = TRUE;
 444             for (k = 0; k < f->order; k++)
 445                 dots[f->dots[k]-g->dots] = TRUE;
 446         }
 447     }
 448
 449     /*
 450      * Work out the new indices of those faces and dots, when we
 451      * compact the arrays containing them.
 452      */
 453     for (i = newfaces = 0; i < g->num_faces; i++)
 454         faces[i] = (faces[i] ? newfaces++ : -1);
 455     for (i = newdots = 0; i < g->num_dots; i++)
 456         dots[i] = (dots[i] ? newdots++ : -1);
 457
 458     /*
 459      * Free the dynamically allocated 'dots' pointer lists in faces
 460      * we're going to discard.
 461      */
 462     for (i = 0; i < g->num_faces; i++)
 463         if (faces[i] < 0)
 464             sfree(g->faces[i].dots);
 465
 466     /*
 467      * Go through and compact the arrays.
 468      */
 469     for (i = 0; i < g->num_dots; i++)
 470         if (dots[i] >= 0) {
 471             grid_dot *dnew = g->dots + dots[i], *dold = g->dots + i;
 472             *dnew = *dold;             /* structure copy */
 473         }
 474     for (i = 0; i < g->num_faces; i++)
 475         if (faces[i] >= 0) {
 476             grid_face *fnew = g->faces + faces[i], *fold = g->faces + i;
 477             *fnew = *fold;             /* structure copy */
 478             for (j = 0; j < fnew->order; j++) {
 479                 /*
 480                  * Reindex the dots in this face.
 481                  */
 482                 k = fnew->dots[j] - g->dots;
 483                 fnew->dots[j] = g->dots + dots[k];
 484             }
 485         }
 486     g->num_faces = newfaces;
 487     g->num_dots = newdots;
 488
 489     sfree(dotpairs);
 490     sfree(dsf);
 491     sfree(dots);
 492     sfree(faces);
 493 }
 494
 495 /* Input: grid has its dots and faces initialised:
 496  * - dots have (optionally) x and y coordinates, but no edges or faces
 497  * (pointers are NULL).
 498  * - edges not initialised at all
 499  * - faces initialised and know which dots they have (but no edges yet).  The
 500  * dots around each face are assumed to be clockwise.
 501  *
 502  * Output: grid is complete and valid with all relationships defined.
 503  */
 504 static void grid_make_consistent(grid *g)
 505 {
 506     int i;
 507     tree234 *incomplete_edges;
 508     grid_edge *next_new_edge; /* Where new edge will go into g->edges */
 509
 510     grid_debug_basic(g);
 511
 512     /* ====== Stage 1 ======
 513      * Generate edges
 514      */
 515
 516     /* We know how many dots and faces there are, so we can find the exact
 517      * number of edges from Euler's polyhedral formula: F + V = E + 2 .
 518      * We use "-1", not "-2" here, because Euler's formula includes the
 519      * infinite face, which we don't count. */
 520     g->num_edges = g->num_faces + g->num_dots - 1;
 521     g->edges = snewn(g->num_edges, grid_edge);
 522     next_new_edge = g->edges;
 523
 524     /* Iterate over faces, and over each face's dots, generating edges as we
 525      * go.  As we find each new edge, we can immediately fill in the edge's
 526      * dots, but only one of the edge's faces.  Later on in the iteration, we
 527      * will find the same edge again (unless it's on the border), but we will
 528      * know the other face.
 529      * For efficiency, maintain a list of the incomplete edges, sorted by
 530      * their dots. */
 531     incomplete_edges = newtree234(grid_edge_bydots_cmpfn);
 532     for (i = 0; i < g->num_faces; i++) {
 533         grid_face *f = g->faces + i;
 534         int j;
 535         for (j = 0; j < f->order; j++) {
 536             grid_edge e; /* fake edge for searching */
 537             grid_edge *edge_found;
 538             int j2 = j + 1;
 539             if (j2 == f->order)
 540                 j2 = 0;
 541             e.dot1 = f->dots[j];
 542             e.dot2 = f->dots[j2];
 543             /* Use del234 instead of find234, because we always want to
 544              * remove the edge if found */
 545             edge_found = del234(incomplete_edges, &e);
 546             if (edge_found) {
 547                 /* This edge already added, so fill out missing face.
 548                  * Edge is already removed from incomplete_edges. */
 549                 edge_found->face2 = f;
 550             } else {
 551                 assert(next_new_edge - g->edges < g->num_edges);
 552                 next_new_edge->dot1 = e.dot1;
 553                 next_new_edge->dot2 = e.dot2;
 554                 next_new_edge->face1 = f;
 555                 next_new_edge->face2 = NULL; /* potentially infinite face */
 556                 add234(incomplete_edges, next_new_edge);
 557                 ++next_new_edge;
 558             }
 559         }
 560     }
 561     freetree234(incomplete_edges);
 562
 563     /* ====== Stage 2 ======
 564      * For each face, build its edge list.
 565      */
 566
 567     /* Allocate space for each edge list.  Can do this, because each face's
 568      * edge-list is the same size as its dot-list. */
 569     for (i = 0; i < g->num_faces; i++) {
 570         grid_face *f = g->faces + i;
 571         int j;
 572         f->edges = snewn(f->order, grid_edge*);
 573         /* Preload with NULLs, to help detect potential bugs. */
 574         for (j = 0; j < f->order; j++)
 575             f->edges[j] = NULL;
 576     }
 577
 578     /* Iterate over each edge, and over both its faces.  Add this edge to
 579      * the face's edge-list, after finding where it should go in the
 580      * sequence. */
 581     for (i = 0; i < g->num_edges; i++) {
 582         grid_edge *e = g->edges + i;
 583         int j;
 584         for (j = 0; j < 2; j++) {
 585             grid_face *f = j ? e->face2 : e->face1;
 586             int k, k2;
 587             if (f == NULL) continue;
 588             /* Find one of the dots around the face */
 589             for (k = 0; k < f->order; k++) {
 590                 if (f->dots[k] == e->dot1)
 591                     break; /* found dot1 */
 592             }
 593             assert(k != f->order); /* Must find the dot around this face */
 594
 595             /* Labelling scheme: as we walk clockwise around the face,
 596              * starting at dot0 (f->dots[0]), we hit:
 597              * (dot0), edge0, dot1, edge1, dot2,...
 598              *
 599              *     0
 600              *  0-----1
 601              *        |
 602              *        |1
 603              *        |
 604              *  3-----2
 605              *     2
 606              *
 607              * Therefore, edgeK joins dotK and dot{K+1}
 608              */
 609
 610             /* Around this face, either the next dot or the previous dot
 611              * must be e->dot2.  Otherwise the edge is wrong. */
 612             k2 = k + 1;
 613             if (k2 == f->order)
 614                 k2 = 0;
 615             if (f->dots[k2] == e->dot2) {
 616                 /* dot1(k) and dot2(k2) go clockwise around this face, so add
 617                  * this edge at position k (see diagram). */
 618                 assert(f->edges[k] == NULL);
 619                 f->edges[k] = e;
 620                 continue;
 621             }
 622             /* Try previous dot */
 623             k2 = k - 1;
 624             if (k2 == -1)
 625                 k2 = f->order - 1;
 626             if (f->dots[k2] == e->dot2) {
 627                 /* dot1(k) and dot2(k2) go anticlockwise around this face. */
 628                 assert(f->edges[k2] == NULL);
 629                 f->edges[k2] = e;
 630                 continue;
 631             }
 632             assert(!"Grid broken: bad edge-face relationship");
 633         }
 634     }
 635
 636     /* ====== Stage 3 ======
 637      * For each dot, build its edge-list and face-list.
 638      */
 639
 640     /* We don't know how many edges/faces go around each dot, so we can't
 641      * allocate the right space for these lists.  Pre-compute the sizes by
 642      * iterating over each edge and recording a tally against each dot. */
 643     for (i = 0; i < g->num_dots; i++) {
 644         g->dots[i].order = 0;
 645     }
 646     for (i = 0; i < g->num_edges; i++) {
 647         grid_edge *e = g->edges + i;
 648         ++(e->dot1->order);
 649         ++(e->dot2->order);
 650     }
 651     /* Now we have the sizes, pre-allocate the edge and face lists. */
 652     for (i = 0; i < g->num_dots; i++) {
 653         grid_dot *d = g->dots + i;
 654         int j;
 655         assert(d->order >= 2); /* sanity check */
 656         d->edges = snewn(d->order, grid_edge*);
 657         d->faces = snewn(d->order, grid_face*);
 658         for (j = 0; j < d->order; j++) {
 659             d->edges[j] = NULL;
 660             d->faces[j] = NULL;
 661         }
 662     }
 663     /* For each dot, need to find a face that touches it, so we can seed
 664      * the edge-face-edge-face process around each dot. */
 665     for (i = 0; i < g->num_faces; i++) {
 666         grid_face *f = g->faces + i;
 667         int j;
 668         for (j = 0; j < f->order; j++) {
 669             grid_dot *d = f->dots[j];
 670             d->faces[0] = f;
 671         }
 672     }
 673     /* Each dot now has a face in its first slot.  Generate the remaining
 674      * faces and edges around the dot, by searching both clockwise and
 675      * anticlockwise from the first face.  Need to do both directions,
 676      * because of the possibility of hitting the infinite face, which
 677      * blocks progress.  But there's only one such face, so we will
 678      * succeed in finding every edge and face this way. */
 679     for (i = 0; i < g->num_dots; i++) {
 680         grid_dot *d = g->dots + i;
 681         int current_face1 = 0; /* ascends clockwise */
 682         int current_face2 = 0; /* descends anticlockwise */
 683
 684         /* Labelling scheme: as we walk clockwise around the dot, starting
 685          * at face0 (d->faces[0]), we hit:
 686          * (face0), edge0, face1, edge1, face2,...
 687          *
 688          *       0
 689          *       |
 690          *    0  |  1
 691          *       |
 692          *  -----d-----1
 693          *       |
 694          *       |  2
 695          *       |
 696          *       2
 697          *
 698          * So, for example, face1 should be joined to edge0 and edge1,
 699          * and those edges should appear in an anticlockwise sense around
 700          * that face (see diagram). */
 701
 702         /* clockwise search */
 703         while (TRUE) {
 704             grid_face *f = d->faces[current_face1];
 705             grid_edge *e;
 706             int j;
 707             assert(f != NULL);
 708             /* find dot around this face */
 709             for (j = 0; j < f->order; j++) {
 710                 if (f->dots[j] == d)
 711                     break;
 712             }
 713             assert(j != f->order); /* must find dot */
 714
 715             /* Around f, required edge is anticlockwise from the dot.  See
 716              * the other labelling scheme higher up, for why we subtract 1
 717              * from j. */
 718             j--;
 719             if (j == -1)
 720                 j = f->order - 1;
 721             e = f->edges[j];
 722             d->edges[current_face1] = e; /* set edge */
 723             current_face1++;
 724             if (current_face1 == d->order)
 725                 break;
 726             else {
 727                 /* set face */
 728                 d->faces[current_face1] =
 729                     (e->face1 == f) ? e->face2 : e->face1;
 730                 if (d->faces[current_face1] == NULL)
 731                     break; /* cannot progress beyond infinite face */
 732             }
 733         }
 734         /* If the clockwise search made it all the way round, don't need to
 735          * bother with the anticlockwise search. */
 736         if (current_face1 == d->order)
 737             continue; /* this dot is complete, move on to next dot */
 738
 739         /* anticlockwise search */
 740         while (TRUE) {
 741             grid_face *f = d->faces[current_face2];
 742             grid_edge *e;
 743             int j;
 744             assert(f != NULL);
 745             /* find dot around this face */
 746             for (j = 0; j < f->order; j++) {
 747                 if (f->dots[j] == d)
 748                     break;
 749             }
 750             assert(j != f->order); /* must find dot */
 751
 752             /* Around f, required edge is clockwise from the dot. */
 753             e = f->edges[j];
 754
 755             current_face2--;
 756             if (current_face2 == -1)
 757                 current_face2 = d->order - 1;
 758             d->edges[current_face2] = e; /* set edge */
 759
 760             /* set face */
 761             if (current_face2 == current_face1)
 762                 break;
 763             d->faces[current_face2] =
 764                     (e->face1 == f) ? e->face2 : e->face1;
 765             /* There's only 1 infinite face, so we must get all the way
 766              * to current_face1 before we hit it. */
 767             assert(d->faces[current_face2]);
 768         }
 769     }
 770
 771     /* ====== Stage 4 ======
 772      * Compute other grid settings
 773      */
 774
 775     /* Bounding rectangle */
 776     for (i = 0; i < g->num_dots; i++) {
 777         grid_dot *d = g->dots + i;
 778         if (i == 0) {
 779             g->lowest_x = g->highest_x = d->x;
 780             g->lowest_y = g->highest_y = d->y;
 781         } else {
 782             g->lowest_x = min(g->lowest_x, d->x);
 783             g->highest_x = max(g->highest_x, d->x);
 784             g->lowest_y = min(g->lowest_y, d->y);
 785             g->highest_y = max(g->highest_y, d->y);
 786         }
 787     }
 788
 789     grid_debug_derived(g);
 790 }
 791
 792 /* Helpers for making grid-generation easier.  These functions are only
 793  * intended for use during grid generation. */
 794
 795 /* Comparison function for the (tree234) sorted dot list */
 796 static int grid_point_cmp_fn(void *v1, void *v2)
 797 {
 798     grid_dot *p1 = v1;
 799     grid_dot *p2 = v2;
 800     if (p1->y != p2->y)
 801         return p2->y - p1->y;
 802     else
 803         return p2->x - p1->x;
 804 }
 805 /* Add a new face to the grid, with its dot list allocated.
 806  * Assumes there's enough space allocated for the new face in grid->faces */
 807 static void grid_face_add_new(grid *g, int face_size)
 808 {
 809     int i;
 810     grid_face *new_face = g->faces + g->num_faces;
 811     new_face->order = face_size;
 812     new_face->dots = snewn(face_size, grid_dot*);
 813     for (i = 0; i < face_size; i++)
 814         new_face->dots[i] = NULL;
 815     new_face->edges = NULL;
 816     new_face->has_incentre = FALSE;
 817     g->num_faces++;
 818 }
 819 /* Assumes dot list has enough space */
 820 static grid_dot *grid_dot_add_new(grid *g, int x, int y)
 821 {
 822     grid_dot *new_dot = g->dots + g->num_dots;
 823     new_dot->order = 0;
 824     new_dot->edges = NULL;
 825     new_dot->faces = NULL;
 826     new_dot->x = x;
 827     new_dot->y = y;
 828     g->num_dots++;
 829     return new_dot;
 830 }
 831 /* Retrieve a dot with these (x,y) coordinates.  Either return an existing dot
 832  * in the dot_list, or add a new dot to the grid (and the dot_list) and
 833  * return that.
 834  * Assumes g->dots has enough capacity allocated */
 835 static grid_dot *grid_get_dot(grid *g, tree234 *dot_list, int x, int y)
 836 {
 837     grid_dot test, *ret;
 838
 839     test.order = 0;
 840     test.edges = NULL;
 841     test.faces = NULL;
 842     test.x = x;
 843     test.y = y;
 844     ret = find234(dot_list, &test, NULL);
 845     if (ret)
 846         return ret;
 847
 848     ret = grid_dot_add_new(g, x, y);
 849     add234(dot_list, ret);
 850     return ret;
 851 }
 852
 853 /* Sets the last face of the grid to include this dot, at this position
 854  * around the face. Assumes num_faces is at least 1 (a new face has
 855  * previously been added, with the required number of dots allocated) */
 856 static void grid_face_set_dot(grid *g, grid_dot *d, int position)
 857 {
 858     grid_face *last_face = g->faces + g->num_faces - 1;
 859     last_face->dots[position] = d;
 860 }
 861
 862 /*
 863  * Helper routines for grid_find_incentre.
 864  */
 865 static int solve_2x2_matrix(double mx[4], double vin[2], double vout[2])
 866 {
 867     double inv[4];
 868     double det;
 869     det = (mx[0]*mx[3] - mx[1]*mx[2]);
 870     if (det == 0)
 871         return FALSE;
 872
 873     inv[0] = mx[3] / det;
 874     inv[1] = -mx[1] / det;
 875     inv[2] = -mx[2] / det;
 876     inv[3] = mx[0] / det;
 877
 878     vout[0] = inv[0]*vin[0] + inv[1]*vin[1];
 879     vout[1] = inv[2]*vin[0] + inv[3]*vin[1];
 880
 881     return TRUE;
 882 }
 883 static int solve_3x3_matrix(double mx[9], double vin[3], double vout[3])
 884 {
 885     double inv[9];
 886     double det;
 887
 888     det = (mx[0]*mx[4]*mx[8] + mx[1]*mx[5]*mx[6] + mx[2]*mx[3]*mx[7] -
 889            mx[0]*mx[5]*mx[7] - mx[1]*mx[3]*mx[8] - mx[2]*mx[4]*mx[6]);
 890     if (det == 0)
 891         return FALSE;
 892
 893     inv[0] = (mx[4]*mx[8] - mx[5]*mx[7]) / det;
 894     inv[1] = (mx[2]*mx[7] - mx[1]*mx[8]) / det;
 895     inv[2] = (mx[1]*mx[5] - mx[2]*mx[4]) / det;
 896     inv[3] = (mx[5]*mx[6] - mx[3]*mx[8]) / det;
 897     inv[4] = (mx[0]*mx[8] - mx[2]*mx[6]) / det;
 898     inv[5] = (mx[2]*mx[3] - mx[0]*mx[5]) / det;
 899     inv[6] = (mx[3]*mx[7] - mx[4]*mx[6]) / det;
 900     inv[7] = (mx[1]*mx[6] - mx[0]*mx[7]) / det;
 901     inv[8] = (mx[0]*mx[4] - mx[1]*mx[3]) / det;
 902
 903     vout[0] = inv[0]*vin[0] + inv[1]*vin[1] + inv[2]*vin[2];
 904     vout[1] = inv[3]*vin[0] + inv[4]*vin[1] + inv[5]*vin[2];
 905     vout[2] = inv[6]*vin[0] + inv[7]*vin[1] + inv[8]*vin[2];
 906
 907     return TRUE;
 908 }
 909
 910 void grid_find_incentre(grid_face *f)
 911 {
 912     double xbest, ybest, bestdist;
 913     int i, j, k, m;
 914     grid_dot *edgedot1[3], *edgedot2[3];
 915     grid_dot *dots[3];
 916     int nedges, ndots;
 917
 918     if (f->has_incentre)
 919         return;
 920
 921     /*
 922      * Find the point in the polygon with the maximum distance to any
 923      * edge or corner.
 924      *
 925      * Such a point must exist which is in contact with at least three
 926      * edges and/or vertices. (Proof: if it's only in contact with two
 927      * edges and/or vertices, it can't even be at a _local_ maximum -
 928      * any such circle can always be expanded in some direction.) So
 929      * we iterate through all 3-subsets of the combined set of edges
 930      * and vertices; for each subset we generate one or two candidate
 931      * points that might be the incentre, and then we vet each one to
 932      * see if it's inside the polygon and what its maximum radius is.
 933      *
 934      * (There's one case which this algorithm will get noticeably
 935      * wrong, and that's when a continuum of equally good answers
 936      * exists due to parallel edges. Consider a long thin rectangle,
 937      * for instance, or a parallelogram. This algorithm will pick a
 938      * point near one end, and choose the end arbitrarily; obviously a
 939      * nicer point to choose would be in the centre. To fix this I
 940      * would have to introduce a special-case system which detected
 941      * parallel edges in advance, set aside all candidate points
 942      * generated using both edges in a parallel pair, and generated
 943      * some additional candidate points half way between them. Also,
 944      * of course, I'd have to cope with rounding error making such a
 945      * point look worse than one of its endpoints. So I haven't done
 946      * this for the moment, and will cross it if necessary when I come
 947      * to it.)
 948      *
 949      * We don't actually iterate literally over _edges_, in the sense
 950      * of grid_edge structures. Instead, we fill in edgedot1[] and
 951      * edgedot2[] with a pair of dots adjacent in the face's list of
 952      * vertices. This ensures that we get the edges in consistent
 953      * orientation, which we could not do from the grid structure
 954      * alone. (A moment's consideration of an order-3 vertex should
 955      * make it clear that if a notional arrow was written on each
 956      * edge, _at least one_ of the three faces bordering that vertex
 957      * would have to have the two arrows tip-to-tip or tail-to-tail
 958      * rather than tip-to-tail.)
 959      */
 960     nedges = ndots = 0;
 961     bestdist = 0;
 962     xbest = ybest = 0;
 963
 964     for (i = 0; i+2 < 2*f->order; i++) {
 965         if (i < f->order) {
 966             edgedot1[nedges] = f->dots[i];
 967             edgedot2[nedges++] = f->dots[(i+1)%f->order];
 968         } else
 969             dots[ndots++] = f->dots[i - f->order];
 970
 971         for (j = i+1; j+1 < 2*f->order; j++) {
 972             if (j < f->order) {
 973                 edgedot1[nedges] = f->dots[j];
 974                 edgedot2[nedges++] = f->dots[(j+1)%f->order];
 975             } else
 976                 dots[ndots++] = f->dots[j - f->order];
 977
 978             for (k = j+1; k < 2*f->order; k++) {
 979                 double cx[2], cy[2];   /* candidate positions */
 980                 int cn = 0;            /* number of candidates */
 981
 982                 if (k < f->order) {
 983                     edgedot1[nedges] = f->dots[k];
 984                     edgedot2[nedges++] = f->dots[(k+1)%f->order];
 985                 } else
 986                     dots[ndots++] = f->dots[k - f->order];
 987
 988                 /*
 989                  * Find a point, or pair of points, equidistant from
 990                  * all the specified edges and/or vertices.
 991                  */
 992                 if (nedges == 3) {
 993                     /*
 994                      * Three edges. This is a linear matrix equation:
 995                      * each row of the matrix represents the fact that
 996                      * the point (x,y) we seek is at distance r from
 997                      * that edge, and we solve three of those
 998                      * simultaneously to obtain x,y,r. (We ignore r.)
 999                      */
1000                     double matrix[9], vector[3], vector2[3];
1001                     int m;
1002
1003                     for (m = 0; m < 3; m++) {
1004                         int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
1005                         int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
1006                         int dx = x2-x1, dy = y2-y1;
1007
1008                         /*
1009                          * ((x,y) - (x1,y1)) . (dy,-dx) = r |(dx,dy)|
1010                          *
1011                          * => x dy - y dx - r |(dx,dy)| = (x1 dy - y1 dx)
1012                          */
1013                         matrix[3*m+0] = dy;
1014                         matrix[3*m+1] = -dx;
1015                         matrix[3*m+2] = -sqrt((double)dx*dx+(double)dy*dy);
1016                         vector[m] = (double)x1*dy - (double)y1*dx;
1017                     }
1018
1019                     if (solve_3x3_matrix(matrix, vector, vector2)) {
1020                         cx[cn] = vector2[0];
1021                         cy[cn] = vector2[1];
1022                         cn++;
1023                     }
1024                 } else if (nedges == 2) {
1025                     /*
1026                      * Two edges and a dot. This will end up in a
1027                      * quadratic equation.
1028                      *
1029                      * First, look at the two edges. Having our point
1030                      * be some distance r from both of them gives rise
1031                      * to a pair of linear equations in x,y,r of the
1032                      * form
1033                      *
1034                      *   (x-x1) dy - (y-y1) dx = r sqrt(dx^2+dy^2)
1035                      *
1036                      * We eliminate r between those equations to give
1037                      * us a single linear equation in x,y describing
1038                      * the locus of points equidistant from both lines
1039                      * - i.e. the angle bisector.
1040                      *
1041                      * We then choose one of x,y to be a parameter t,
1042                      * and derive linear formulae for x,y,r in terms
1043                      * of t. This enables us to write down the
1044                      * circular equation (x-xd)^2+(y-yd)^2=r^2 as a
1045                      * quadratic in t; solving that and substituting
1046                      * in for x,y gives us two candidate points.
1047                      */
1048                     double eqs[2][4];  /* a,b,c,d : ax+by+cr=d */
1049                     double eq[3];      /* a,b,c: ax+by=c */
1050                     double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
1051                     double q[3];                /* a,b,c: at^2+bt+c=0 */
1052                     double disc;
1053
1054                     /* Find equations of the two input lines. */
1055                     for (m = 0; m < 2; m++) {
1056                         int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
1057                         int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
1058                         int dx = x2-x1, dy = y2-y1;
1059
1060                         eqs[m][0] = dy;
1061                         eqs[m][1] = -dx;
1062                         eqs[m][2] = -sqrt(dx*dx+dy*dy);
1063                         eqs[m][3] = x1*dy - y1*dx;
1064                     }
1065
1066                     /* Derive the angle bisector by eliminating r. */
1067                     eq[0] = eqs[0][0]*eqs[1][2] - eqs[1][0]*eqs[0][2];
1068                     eq[1] = eqs[0][1]*eqs[1][2] - eqs[1][1]*eqs[0][2];
1069                     eq[2] = eqs[0][3]*eqs[1][2] - eqs[1][3]*eqs[0][2];
1070
1071                     /* Parametrise x and y in terms of some t. */
1072                     if (fabs(eq[0]) < fabs(eq[1])) {
1073                         /* Parameter is x. */
1074                         xt[0] = 1; xt[1] = 0;
1075                         yt[0] = -eq[0]/eq[1]; yt[1] = eq[2]/eq[1];
1076                     } else {
1077                         /* Parameter is y. */
1078                         yt[0] = 1; yt[1] = 0;
1079                         xt[0] = -eq[1]/eq[0]; xt[1] = eq[2]/eq[0];
1080                     }
1081
1082                     /* Find a linear representation of r using eqs[0]. */
1083                     rt[0] = -(eqs[0][0]*xt[0] + eqs[0][1]*yt[0])/eqs[0][2];
1084                     rt[1] = (eqs[0][3] - eqs[0][0]*xt[1] -
1085                              eqs[0][1]*yt[1])/eqs[0][2];
1086
1087                     /* Construct the quadratic equation. */
1088                     q[0] = -rt[0]*rt[0];
1089                     q[1] = -2*rt[0]*rt[1];
1090                     q[2] = -rt[1]*rt[1];
1091                     q[0] += xt[0]*xt[0];
1092                     q[1] += 2*xt[0]*(xt[1]-dots[0]->x);
1093                     q[2] += (xt[1]-dots[0]->x)*(xt[1]-dots[0]->x);
1094                     q[0] += yt[0]*yt[0];
1095                     q[1] += 2*yt[0]*(yt[1]-dots[0]->y);
1096                     q[2] += (yt[1]-dots[0]->y)*(yt[1]-dots[0]->y);
1097
1098                     /* And solve it. */
1099                     disc = q[1]*q[1] - 4*q[0]*q[2];
1100                     if (disc >= 0) {
1101                         double t;
1102
1103                         disc = sqrt(disc);
1104
1105                         t = (-q[1] + disc) / (2*q[0]);
1106                         cx[cn] = xt[0]*t + xt[1];
1107                         cy[cn] = yt[0]*t + yt[1];
1108                         cn++;
1109
1110                         t = (-q[1] - disc) / (2*q[0]);
1111                         cx[cn] = xt[0]*t + xt[1];
1112                         cy[cn] = yt[0]*t + yt[1];
1113                         cn++;
1114                     }
1115                 } else if (nedges == 1) {
1116                     /*
1117                      * Two dots and an edge. This one's another
1118                      * quadratic equation.
1119                      *
1120                      * The point we want must lie on the perpendicular
1121                      * bisector of the two dots; that much is obvious.
1122                      * So we can construct a parametrisation of that
1123                      * bisecting line, giving linear formulae for x,y
1124                      * in terms of t. We can also express the distance
1125                      * from the edge as such a linear formula.
1126                      *
1127                      * Then we set that equal to the radius of the
1128                      * circle passing through the two points, which is
1129                      * a Pythagoras exercise; that gives rise to a
1130                      * quadratic in t, which we solve.
1131                      */
1132                     double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
1133                     double q[3];                /* a,b,c: at^2+bt+c=0 */
1134                     double disc;
1135                     double halfsep;
1136
1137                     /* Find parametric formulae for x,y. */
1138                     {
1139                         int x1 = dots[0]->x, x2 = dots[1]->x;
1140                         int y1 = dots[0]->y, y2 = dots[1]->y;
1141                         int dx = x2-x1, dy = y2-y1;
1142                         double d = sqrt((double)dx*dx + (double)dy*dy);
1143
1144                         xt[1] = (x1+x2)/2.0;
1145                         yt[1] = (y1+y2)/2.0;
1146                         /* It's convenient if we have t at standard scale. */
1147                         xt[0] = -dy/d;
1148                         yt[0] = dx/d;
1149
1150                         /* Also note down half the separation between
1151                          * the dots, for use in computing the circle radius. */
1152                         halfsep = 0.5*d;
1153                     }
1154
1155                     /* Find a parametric formula for r. */
1156                     {
1157                         int x1 = edgedot1[0]->x, x2 = edgedot2[0]->x;
1158                         int y1 = edgedot1[0]->y, y2 = edgedot2[0]->y;
1159                         int dx = x2-x1, dy = y2-y1;
1160                         double d = sqrt((double)dx*dx + (double)dy*dy);
1161                         rt[0] = (xt[0]*dy - yt[0]*dx) / d;
1162                         rt[1] = ((xt[1]-x1)*dy - (yt[1]-y1)*dx) / d;
1163                     }
1164
1165                     /* Construct the quadratic equation. */
1166                     q[0] = rt[0]*rt[0];
1167                     q[1] = 2*rt[0]*rt[1];
1168                     q[2] = rt[1]*rt[1];
1169                     q[0] -= 1;
1170                     q[2] -= halfsep*halfsep;
1171
1172                     /* And solve it. */
1173                     disc = q[1]*q[1] - 4*q[0]*q[2];
1174                     if (disc >= 0) {
1175                         double t;
1176
1177                         disc = sqrt(disc);
1178
1179                         t = (-q[1] + disc) / (2*q[0]);
1180                         cx[cn] = xt[0]*t + xt[1];
1181                         cy[cn] = yt[0]*t + yt[1];
1182                         cn++;
1183
1184                         t = (-q[1] - disc) / (2*q[0]);
1185                         cx[cn] = xt[0]*t + xt[1];
1186                         cy[cn] = yt[0]*t + yt[1];
1187                         cn++;
1188                     }
1189                 } else if (nedges == 0) {
1190                     /*
1191                      * Three dots. This is another linear matrix
1192                      * equation, this time with each row of the matrix
1193                      * representing the perpendicular bisector between
1194                      * two of the points. Of course we only need two
1195                      * such lines to find their intersection, so we
1196                      * need only solve a 2x2 matrix equation.
1197                      */
1198
1199                     double matrix[4], vector[2], vector2[2];
1200                     int m;
1201
1202                     for (m = 0; m < 2; m++) {
1203                         int x1 = dots[m]->x, x2 = dots[m+1]->x;
1204                         int y1 = dots[m]->y, y2 = dots[m+1]->y;
1205                         int dx = x2-x1, dy = y2-y1;
1206
1207                         /*
1208                          * ((x,y) - (x1,y1)) . (dx,dy) = 1/2 |(dx,dy)|^2
1209                          *
1210                          * => 2x dx + 2y dy = dx^2+dy^2 + (2 x1 dx + 2 y1 dy)
1211                          */
1212                         matrix[2*m+0] = 2*dx;
1213                         matrix[2*m+1] = 2*dy;
1214                         vector[m] = ((double)dx*dx + (double)dy*dy +
1215                                      2.0*x1*dx + 2.0*y1*dy);
1216                     }
1217
1218                     if (solve_2x2_matrix(matrix, vector, vector2)) {
1219                         cx[cn] = vector2[0];
1220                         cy[cn] = vector2[1];
1221                         cn++;
1222                     }
1223                 }
1224
1225                 /*
1226                  * Now go through our candidate points and see if any
1227                  * of them are better than what we've got so far.
1228                  */
1229                 for (m = 0; m < cn; m++) {
1230                     double x = cx[m], y = cy[m];
1231
1232                     /*
1233                      * First, disqualify the point if it's not inside
1234                      * the polygon, which we work out by counting the
1235                      * edges to the right of the point. (For
1236                      * tiebreaking purposes when edges start or end on
1237                      * our y-coordinate or go right through it, we
1238                      * consider our point to be offset by a small
1239                      * _positive_ epsilon in both the x- and
1240                      * y-direction.)
1241                      */
1242                     int e, in = 0;
1243                     for (e = 0; e < f->order; e++) {
1244                         int xs = f->edges[e]->dot1->x;
1245                         int xe = f->edges[e]->dot2->x;
1246                         int ys = f->edges[e]->dot1->y;
1247                         int ye = f->edges[e]->dot2->y;
1248                         if ((y >= ys && y < ye) || (y >= ye && y < ys)) {
1249                             /*
1250                              * The line goes past our y-position. Now we need
1251                              * to know if its x-coordinate when it does so is
1252                              * to our right.
1253                              *
1254                              * The x-coordinate in question is mathematically
1255                              * (y - ys) * (xe - xs) / (ye - ys), and we want
1256                              * to know whether (x - xs) >= that. Of course we
1257                              * avoid the division, so we can work in integers;
1258                              * to do this we must multiply both sides of the
1259                              * inequality by ye - ys, which means we must
1260                              * first check that's not negative.
1261                              */
1262                             int num = xe - xs, denom = ye - ys;
1263                             if (denom < 0) {
1264                                 num = -num;
1265                                 denom = -denom;
1266                             }
1267                             if ((x - xs) * denom >= (y - ys) * num)
1268                                 in ^= 1;
1269                         }
1270                     }
1271
1272                     if (in) {
1273 #ifdef HUGE_VAL
1274                         double mindist = HUGE_VAL;
1275 #else
1276 #ifdef DBL_MAX
1277                         double mindist = DBL_MAX;
1278 #else
1279 #error No way to get maximum floating-point number.
1280 #endif
1281 #endif
1282                         int e, d;
1283
1284                         /*
1285                          * This point is inside the polygon, so now we check
1286                          * its minimum distance to every edge and corner.
1287                          * First the corners ...
1288                          */
1289                         for (d = 0; d < f->order; d++) {
1290                             int xp = f->dots[d]->x;
1291                             int yp = f->dots[d]->y;
1292                             double dx = x - xp, dy = y - yp;
1293                             double dist = dx*dx + dy*dy;
1294                             if (mindist > dist)
1295                                 mindist = dist;
1296                         }
1297
1298                         /*
1299                          * ... and now also check the perpendicular distance
1300                          * to every edge, if the perpendicular lies between
1301                          * the edge's endpoints.
1302                          */
1303                         for (e = 0; e < f->order; e++) {
1304                             int xs = f->edges[e]->dot1->x;
1305                             int xe = f->edges[e]->dot2->x;
1306                             int ys = f->edges[e]->dot1->y;
1307                             int ye = f->edges[e]->dot2->y;
1308
1309                             /*
1310                              * If s and e are our endpoints, and p our
1311                              * candidate circle centre, the foot of a
1312                              * perpendicular from p to the line se lies
1313                              * between s and e if and only if (p-s).(e-s) lies
1314                              * strictly between 0 and (e-s).(e-s).
1315                              */
1316                             int edx = xe - xs, edy = ye - ys;
1317                             double pdx = x - xs, pdy = y - ys;
1318                             double pde = pdx * edx + pdy * edy;
1319                             long ede = (long)edx * edx + (long)edy * edy;
1320                             if (0 < pde && pde < ede) {
1321                                 /*
1322                                  * Yes, the nearest point on this edge is
1323                                  * closer than either endpoint, so we must
1324                                  * take it into account by measuring the
1325                                  * perpendicular distance to the edge and
1326                                  * checking its square against mindist.
1327                                  */
1328
1329                                 double pdre = pdx * edy - pdy * edx;
1330                                 double sqlen = pdre * pdre / ede;
1331
1332                                 if (mindist > sqlen)
1333                                     mindist = sqlen;
1334                             }
1335                         }
1336
1337                         /*
1338                          * Right. Now we know the biggest circle around this
1339                          * point, so we can check it against bestdist.
1340                          */
1341                         if (bestdist < mindist) {
1342                             bestdist = mindist;
1343                             xbest = x;
1344                             ybest = y;
1345                         }
1346                     }
1347                 }
1348
1349                 if (k < f->order)
1350                     nedges--;
1351                 else
1352                     ndots--;
1353             }
1354             if (j < f->order)
1355                 nedges--;
1356             else
1357                 ndots--;
1358         }
1359         if (i < f->order)
1360             nedges--;
1361         else
1362             ndots--;
1363     }
1364
1365     assert(bestdist > 0);
1366
1367     f->has_incentre = TRUE;
1368     f->ix = xbest + 0.5;               /* round to nearest */
1369     f->iy = ybest + 0.5;
1370 }
1371
1372 /* ------ Generate various types of grid ------ */
1373
1374 /* General method is to generate faces, by calculating their dot coordinates.
1375  * As new faces are added, we keep track of all the dots so we can tell when
1376  * a new face reuses an existing dot.  For example, two squares touching at an
1377  * edge would generate six unique dots: four dots from the first face, then
1378  * two additional dots for the second face, because we detect the other two
1379  * dots have already been taken up.  This list is stored in a tree234
1380  * called "points".  No extra memory-allocation needed here - we store the
1381  * actual grid_dot* pointers, which all point into the g->dots list.
1382  * For this reason, we have to calculate coordinates in such a way as to
1383  * eliminate any rounding errors, so we can detect when a dot on one
1384  * face precisely lands on a dot of a different face.  No floating-point
1385  * arithmetic here!
1386  */
1387
1388 #define SQUARE_TILESIZE 20
1389
1390 static void grid_size_square(int width, int height,
1391                       int *tilesize, int *xextent, int *yextent)
1392 {
1393     int a = SQUARE_TILESIZE;
1394
1395     *tilesize = a;
1396     *xextent = width * a;
1397     *yextent = height * a;
1398 }
1399
1400 static grid *grid_new_square(int width, int height, const char *desc)
1401 {
1402     int x, y;
1403     /* Side length */
1404     int a = SQUARE_TILESIZE;
1405
1406     /* Upper bounds - don't have to be exact */
1407     int max_faces = width * height;
1408     int max_dots = (width + 1) * (height + 1);
1409
1410     tree234 *points;
1411
1412     grid *g = grid_empty();
1413     g->tilesize = a;
1414     g->faces = snewn(max_faces, grid_face);
1415     g->dots = snewn(max_dots, grid_dot);
1416
1417     points = newtree234(grid_point_cmp_fn);
1418
1419     /* generate square faces */
1420     for (y = 0; y < height; y++) {
1421         for (x = 0; x < width; x++) {
1422             grid_dot *d;
1423             /* face position */
1424             int px = a * x;
1425             int py = a * y;
1426
1427             grid_face_add_new(g, 4);
1428             d = grid_get_dot(g, points, px, py);
1429             grid_face_set_dot(g, d, 0);
1430             d = grid_get_dot(g, points, px + a, py);
1431             grid_face_set_dot(g, d, 1);
1432             d = grid_get_dot(g, points, px + a, py + a);
1433             grid_face_set_dot(g, d, 2);
1434             d = grid_get_dot(g, points, px, py + a);
1435             grid_face_set_dot(g, d, 3);
1436         }
1437     }
1438
1439     freetree234(points);
1440     assert(g->num_faces <= max_faces);
1441     assert(g->num_dots <= max_dots);
1442
1443     grid_make_consistent(g);
1444     return g;
1445 }
1446
1447 #define HONEY_TILESIZE 45
1448 /* Vector for side of hexagon - ratio is close to sqrt(3) */
1449 #define HONEY_A 15
1450 #define HONEY_B 26
1451
1452 static void grid_size_honeycomb(int width, int height,
1453                          int *tilesize, int *xextent, int *yextent)
1454 {
1455     int a = HONEY_A;
1456     int b = HONEY_B;
1457
1458     *tilesize = HONEY_TILESIZE;
1459     *xextent = (3 * a * (width-1)) + 4*a;
1460     *yextent = (2 * b * (height-1)) + 3*b;
1461 }
1462
1463 static grid *grid_new_honeycomb(int width, int height, const char *desc)
1464 {
1465     int x, y;
1466     int a = HONEY_A;
1467     int b = HONEY_B;
1468
1469     /* Upper bounds - don't have to be exact */
1470     int max_faces = width * height;
1471     int max_dots = 2 * (width + 1) * (height + 1);
1472
1473     tree234 *points;
1474
1475     grid *g = grid_empty();
1476     g->tilesize = HONEY_TILESIZE;
1477     g->faces = snewn(max_faces, grid_face);
1478     g->dots = snewn(max_dots, grid_dot);
1479
1480     points = newtree234(grid_point_cmp_fn);
1481
1482     /* generate hexagonal faces */
1483     for (y = 0; y < height; y++) {
1484         for (x = 0; x < width; x++) {
1485             grid_dot *d;
1486             /* face centre */
1487             int cx = 3 * a * x;
1488             int cy = 2 * b * y;
1489             if (x % 2)
1490                 cy += b;
1491             grid_face_add_new(g, 6);
1492
1493             d = grid_get_dot(g, points, cx - a, cy - b);
1494             grid_face_set_dot(g, d, 0);
1495             d = grid_get_dot(g, points, cx + a, cy - b);
1496             grid_face_set_dot(g, d, 1);
1497             d = grid_get_dot(g, points, cx + 2*a, cy);
1498             grid_face_set_dot(g, d, 2);
1499             d = grid_get_dot(g, points, cx + a, cy + b);
1500             grid_face_set_dot(g, d, 3);
1501             d = grid_get_dot(g, points, cx - a, cy + b);
1502             grid_face_set_dot(g, d, 4);
1503             d = grid_get_dot(g, points, cx - 2*a, cy);
1504             grid_face_set_dot(g, d, 5);
1505         }
1506     }
1507
1508     freetree234(points);
1509     assert(g->num_faces <= max_faces);
1510     assert(g->num_dots <= max_dots);
1511
1512     grid_make_consistent(g);
1513     return g;
1514 }
1515
1516 #define TRIANGLE_TILESIZE 18
1517 /* Vector for side of triangle - ratio is close to sqrt(3) */
1518 #define TRIANGLE_VEC_X 15
1519 #define TRIANGLE_VEC_Y 26
1520
1521 static void grid_size_triangular(int width, int height,
1522                           int *tilesize, int *xextent, int *yextent)
1523 {
1524     int vec_x = TRIANGLE_VEC_X;
1525     int vec_y = TRIANGLE_VEC_Y;
1526
1527     *tilesize = TRIANGLE_TILESIZE;
1528     *xextent = (width+1) * 2 * vec_x;
1529     *yextent = height * vec_y;
1530 }
1531
1532 static char *grid_validate_desc_triangular(grid_type type, int width,
1533                                            int height, const char *desc)
1534 {
1535     /*
1536      * Triangular grids: an absent description is valid (indicating
1537      * the old-style approach which had 'ears', i.e. triangles
1538      * connected to only one other face, at some grid corners), and so
1539      * is a description reading just "0" (indicating the new-style
1540      * approach in which those ears are trimmed off). Anything else is
1541      * illegal.
1542      */
1543
1544     if (!desc || !strcmp(desc, "0"))
1545         return NULL;
1546
1547     return "Unrecognised grid description.";
1548 }
1549
1550 /* Doesn't use the previous method of generation, it pre-dates it!
1551  * A triangular grid is just about simple enough to do by "brute force" */
1552 static grid *grid_new_triangular(int width, int height, const char *desc)
1553 {
1554     int x,y;
1555     int version = (desc == NULL ? -1 : atoi(desc));
1556
1557     /* Vector for side of triangle - ratio is close to sqrt(3) */
1558     int vec_x = TRIANGLE_VEC_X;
1559     int vec_y = TRIANGLE_VEC_Y;
1560
1561     int index;
1562
1563     /* convenient alias */
1564     int w = width + 1;
1565
1566     grid *g = grid_empty();
1567     g->tilesize = TRIANGLE_TILESIZE;
1568
1569     if (version == -1) {
1570         /*
1571          * Old-style triangular grid generation, preserved as-is for
1572          * backwards compatibility with old game ids, in which it's
1573          * just a little asymmetric and there are 'ears' (faces linked
1574          * to only one other face) at two grid corners.
1575          *
1576          * Example old-style game ids, which should still work, and in
1577          * which you should see the ears in the TL/BR corners on the
1578          * first one and in the TL/BL corners on the second:
1579          *
1580          *   5x5t1:2c2a1a2a201a1a1a1112a1a2b1211f0b21a2a2a0a
1581          *   5x6t1:a022a212h1a1d1a12c2b11a012b1a20d1a0a12e
1582          */
1583
1584         g->num_faces = width * height * 2;
1585         g->num_dots = (width + 1) * (height + 1);
1586         g->faces = snewn(g->num_faces, grid_face);
1587         g->dots = snewn(g->num_dots, grid_dot);
1588
1589         /* generate dots */
1590         index = 0;
1591         for (y = 0; y <= height; y++) {
1592             for (x = 0; x <= width; x++) {
1593                 grid_dot *d = g->dots + index;
1594                 /* odd rows are offset to the right */
1595                 d->order = 0;
1596                 d->edges = NULL;
1597                 d->faces = NULL;
1598                 d->x = x * 2 * vec_x + ((y % 2) ? vec_x : 0);
1599                 d->y = y * vec_y;
1600                 index++;
1601             }
1602         }
1603
1604         /* generate faces */
1605         index = 0;
1606         for (y = 0; y < height; y++) {
1607             for (x = 0; x < width; x++) {
1608                 /* initialise two faces for this (x,y) */
1609                 grid_face *f1 = g->faces + index;
1610                 grid_face *f2 = f1 + 1;
1611                 f1->edges = NULL;
1612                 f1->order = 3;
1613                 f1->dots = snewn(f1->order, grid_dot*);
1614                 f1->has_incentre = FALSE;
1615                 f2->edges = NULL;
1616                 f2->order = 3;
1617                 f2->dots = snewn(f2->order, grid_dot*);
1618                 f2->has_incentre = FALSE;
1619
1620                 /* face descriptions depend on whether the row-number is
1621                  * odd or even */
1622                 if (y % 2) {
1623                     f1->dots[0] = g->dots + y       * w + x;
1624                     f1->dots[1] = g->dots + (y + 1) * w + x + 1;
1625                     f1->dots[2] = g->dots + (y + 1) * w + x;
1626                     f2->dots[0] = g->dots + y       * w + x;
1627                     f2->dots[1] = g->dots + y       * w + x + 1;
1628                     f2->dots[2] = g->dots + (y + 1) * w + x + 1;
1629                 } else {
1630                     f1->dots[0] = g->dots + y       * w + x;
1631                     f1->dots[1] = g->dots + y       * w + x + 1;
1632                     f1->dots[2] = g->dots + (y + 1) * w + x;
1633                     f2->dots[0] = g->dots + y       * w + x + 1;
1634                     f2->dots[1] = g->dots + (y + 1) * w + x + 1;
1635                     f2->dots[2] = g->dots + (y + 1) * w + x;
1636                 }
1637                 index += 2;
1638             }
1639         }
1640     } else {
1641         /*
1642          * New-style approach, in which there are never any 'ears',
1643          * and if height is even then the grid is nicely 4-way
1644          * symmetric.
1645          *
1646          * Example new-style grids:
1647          *
1648          *   5x5t1:0_21120b11a1a01a1a00c1a0b211021c1h1a2a1a0a
1649          *   5x6t1:0_a1212c22c2a02a2f22a0c12a110d0e1c0c0a101121a1
1650          */
1651         tree234 *points = newtree234(grid_point_cmp_fn);
1652         /* Upper bounds - don't have to be exact */
1653         int max_faces = height * (2*width+1);
1654         int max_dots = (height+1) * (width+1) * 4;
1655
1656         g->faces = snewn(max_faces, grid_face);
1657         g->dots = snewn(max_dots, grid_dot);
1658
1659         for (y = 0; y < height; y++) {
1660             /*
1661              * Each row contains (width+1) triangles one way up, and
1662              * (width) triangles the other way up. Which way up is
1663              * which varies with parity of y. Also, the dots around
1664              * each face will flip direction with parity of y, so we
1665              * set up n1 and n2 to cope with that easily.
1666              */
1667             int y0, y1, n1, n2;
1668             y0 = y1 = y * vec_y;
1669             if (y % 2) {
1670                 y1 += vec_y;
1671                 n1 = 2; n2 = 1;
1672             } else {
1673                 y0 += vec_y;
1674                 n1 = 1; n2 = 2;
1675             }
1676
1677             for (x = 0; x <= width; x++) {
1678                 int x0 = 2*x * vec_x, x1 = x0 + vec_x, x2 = x1 + vec_x;
1679
1680                 /*
1681                  * If the grid has odd height, then we skip the first
1682                  * and last triangles on this row, otherwise they'll
1683                  * end up as ears.
1684                  */
1685                 if (height % 2 == 1 && y == height-1 && (x == 0 || x == width))
1686                     continue;
1687
1688                 grid_face_add_new(g, 3);
1689                 grid_face_set_dot(g, grid_get_dot(g, points, x0, y0), 0);
1690                 grid_face_set_dot(g, grid_get_dot(g, points, x1, y1), n1);
1691                 grid_face_set_dot(g, grid_get_dot(g, points, x2, y0), n2);
1692             }
1693
1694             for (x = 0; x < width; x++) {
1695                 int x0 = (2*x+1) * vec_x, x1 = x0 + vec_x, x2 = x1 + vec_x;
1696
1697                 grid_face_add_new(g, 3);
1698                 grid_face_set_dot(g, grid_get_dot(g, points, x0, y1), 0);
1699                 grid_face_set_dot(g, grid_get_dot(g, points, x1, y0), n2);
1700                 grid_face_set_dot(g, grid_get_dot(g, points, x2, y1), n1);
1701             }
1702         }
1703
1704         freetree234(points);
1705         assert(g->num_faces <= max_faces);
1706         assert(g->num_dots <= max_dots);
1707     }
1708
1709     grid_make_consistent(g);
1710     return g;
1711 }
1712
1713 #define SNUBSQUARE_TILESIZE 18
1714 /* Vector for side of triangle - ratio is close to sqrt(3) */
1715 #define SNUBSQUARE_A 15
1716 #define SNUBSQUARE_B 26
1717
1718 static void grid_size_snubsquare(int width, int height,
1719                           int *tilesize, int *xextent, int *yextent)
1720 {
1721     int a = SNUBSQUARE_A;
1722     int b = SNUBSQUARE_B;
1723
1724     *tilesize = SNUBSQUARE_TILESIZE;
1725     *xextent = (a+b) * (width-1) + a + b;
1726     *yextent = (a+b) * (height-1) + a + b;
1727 }
1728
1729 static grid *grid_new_snubsquare(int width, int height, const char *desc)
1730 {
1731     int x, y;
1732     int a = SNUBSQUARE_A;
1733     int b = SNUBSQUARE_B;
1734
1735     /* Upper bounds - don't have to be exact */
1736     int max_faces = 3 * width * height;
1737     int max_dots = 2 * (width + 1) * (height + 1);
1738
1739     tree234 *points;
1740
1741     grid *g = grid_empty();
1742     g->tilesize = SNUBSQUARE_TILESIZE;
1743     g->faces = snewn(max_faces, grid_face);
1744     g->dots = snewn(max_dots, grid_dot);
1745
1746     points = newtree234(grid_point_cmp_fn);
1747
1748     for (y = 0; y < height; y++) {
1749         for (x = 0; x < width; x++) {
1750             grid_dot *d;
1751             /* face position */
1752             int px = (a + b) * x;
1753             int py = (a + b) * y;
1754
1755             /* generate square faces */
1756             grid_face_add_new(g, 4);
1757             if ((x + y) % 2) {
1758                 d = grid_get_dot(g, points, px + a, py);
1759                 grid_face_set_dot(g, d, 0);
1760                 d = grid_get_dot(g, points, px + a + b, py + a);
1761                 grid_face_set_dot(g, d, 1);
1762                 d = grid_get_dot(g, points, px + b, py + a + b);
1763                 grid_face_set_dot(g, d, 2);
1764                 d = grid_get_dot(g, points, px, py + b);
1765                 grid_face_set_dot(g, d, 3);
1766             } else {
1767                 d = grid_get_dot(g, points, px + b, py);
1768                 grid_face_set_dot(g, d, 0);
1769                 d = grid_get_dot(g, points, px + a + b, py + b);
1770                 grid_face_set_dot(g, d, 1);
1771                 d = grid_get_dot(g, points, px + a, py + a + b);
1772                 grid_face_set_dot(g, d, 2);
1773                 d = grid_get_dot(g, points, px, py + a);
1774                 grid_face_set_dot(g, d, 3);
1775             }
1776
1777             /* generate up/down triangles */
1778             if (x > 0) {
1779                 grid_face_add_new(g, 3);
1780                 if ((x + y) % 2) {
1781                     d = grid_get_dot(g, points, px + a, py);
1782                     grid_face_set_dot(g, d, 0);
1783                     d = grid_get_dot(g, points, px, py + b);
1784                     grid_face_set_dot(g, d, 1);
1785                     d = grid_get_dot(g, points, px - a, py);
1786                     grid_face_set_dot(g, d, 2);
1787                 } else {
1788                     d = grid_get_dot(g, points, px, py + a);
1789                     grid_face_set_dot(g, d, 0);
1790                     d = grid_get_dot(g, points, px + a, py + a + b);
1791                     grid_face_set_dot(g, d, 1);
1792                     d = grid_get_dot(g, points, px - a, py + a + b);
1793                     grid_face_set_dot(g, d, 2);
1794                 }
1795             }
1796
1797             /* generate left/right triangles */
1798             if (y > 0) {
1799                 grid_face_add_new(g, 3);
1800                 if ((x + y) % 2) {
1801                     d = grid_get_dot(g, points, px + a, py);
1802                     grid_face_set_dot(g, d, 0);
1803                     d = grid_get_dot(g, points, px + a + b, py - a);
1804                     grid_face_set_dot(g, d, 1);
1805                     d = grid_get_dot(g, points, px + a + b, py + a);
1806                     grid_face_set_dot(g, d, 2);
1807                 } else {
1808                     d = grid_get_dot(g, points, px, py - a);
1809                     grid_face_set_dot(g, d, 0);
1810                     d = grid_get_dot(g, points, px + b, py);
1811                     grid_face_set_dot(g, d, 1);
1812                     d = grid_get_dot(g, points, px, py + a);
1813                     grid_face_set_dot(g, d, 2);
1814                 }
1815             }
1816         }
1817     }
1818
1819     freetree234(points);
1820     assert(g->num_faces <= max_faces);
1821     assert(g->num_dots <= max_dots);
1822
1823     grid_make_consistent(g);
1824     return g;
1825 }
1826
1827 #define CAIRO_TILESIZE 40
1828 /* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
1829 #define CAIRO_A 14
1830 #define CAIRO_B 31
1831
1832 static void grid_size_cairo(int width, int height,
1833                           int *tilesize, int *xextent, int *yextent)
1834 {
1835     int b = CAIRO_B; /* a unused in determining grid size. */
1836
1837     *tilesize = CAIRO_TILESIZE;
1838     *xextent = 2*b*(width-1) + 2*b;
1839     *yextent = 2*b*(height-1) + 2*b;
1840 }
1841
1842 static grid *grid_new_cairo(int width, int height, const char *desc)
1843 {
1844     int x, y;
1845     int a = CAIRO_A;
1846     int b = CAIRO_B;
1847
1848     /* Upper bounds - don't have to be exact */
1849     int max_faces = 2 * width * height;
1850     int max_dots = 3 * (width + 1) * (height + 1);
1851
1852     tree234 *points;
1853
1854     grid *g = grid_empty();
1855     g->tilesize = CAIRO_TILESIZE;
1856     g->faces = snewn(max_faces, grid_face);
1857     g->dots = snewn(max_dots, grid_dot);
1858
1859     points = newtree234(grid_point_cmp_fn);
1860
1861     for (y = 0; y < height; y++) {
1862         for (x = 0; x < width; x++) {
1863             grid_dot *d;
1864             /* cell position */
1865             int px = 2 * b * x;
1866             int py = 2 * b * y;
1867
1868             /* horizontal pentagons */
1869             if (y > 0) {
1870                 grid_face_add_new(g, 5);
1871                 if ((x + y) % 2) {
1872                     d = grid_get_dot(g, points, px + a, py - b);
1873                     grid_face_set_dot(g, d, 0);
1874                     d = grid_get_dot(g, points, px + 2*b - a, py - b);
1875                     grid_face_set_dot(g, d, 1);
1876                     d = grid_get_dot(g, points, px + 2*b, py);
1877                     grid_face_set_dot(g, d, 2);
1878                     d = grid_get_dot(g, points, px + b, py + a);
1879                     grid_face_set_dot(g, d, 3);
1880                     d = grid_get_dot(g, points, px, py);
1881                     grid_face_set_dot(g, d, 4);
1882                 } else {
1883                     d = grid_get_dot(g, points, px, py);
1884                     grid_face_set_dot(g, d, 0);
1885                     d = grid_get_dot(g, points, px + b, py - a);
1886                     grid_face_set_dot(g, d, 1);
1887                     d = grid_get_dot(g, points, px + 2*b, py);
1888                     grid_face_set_dot(g, d, 2);
1889                     d = grid_get_dot(g, points, px + 2*b - a, py + b);
1890                     grid_face_set_dot(g, d, 3);
1891                     d = grid_get_dot(g, points, px + a, py + b);
1892                     grid_face_set_dot(g, d, 4);
1893                 }
1894             }
1895             /* vertical pentagons */
1896             if (x > 0) {
1897                 grid_face_add_new(g, 5);
1898                 if ((x + y) % 2) {
1899                     d = grid_get_dot(g, points, px, py);
1900                     grid_face_set_dot(g, d, 0);
1901                     d = grid_get_dot(g, points, px + b, py + a);
1902                     grid_face_set_dot(g, d, 1);
1903                     d = grid_get_dot(g, points, px + b, py + 2*b - a);
1904                     grid_face_set_dot(g, d, 2);
1905                     d = grid_get_dot(g, points, px, py + 2*b);
1906                     grid_face_set_dot(g, d, 3);
1907                     d = grid_get_dot(g, points, px - a, py + b);
1908                     grid_face_set_dot(g, d, 4);
1909                 } else {
1910                     d = grid_get_dot(g, points, px, py);
1911                     grid_face_set_dot(g, d, 0);
1912                     d = grid_get_dot(g, points, px + a, py + b);
1913                     grid_face_set_dot(g, d, 1);
1914                     d = grid_get_dot(g, points, px, py + 2*b);
1915                     grid_face_set_dot(g, d, 2);
1916                     d = grid_get_dot(g, points, px - b, py + 2*b - a);
1917                     grid_face_set_dot(g, d, 3);
1918                     d = grid_get_dot(g, points, px - b, py + a);
1919                     grid_face_set_dot(g, d, 4);
1920                 }
1921             }
1922         }
1923     }
1924
1925     freetree234(points);
1926     assert(g->num_faces <= max_faces);
1927     assert(g->num_dots <= max_dots);
1928
1929     grid_make_consistent(g);
1930     return g;
1931 }
1932
1933 #define GREATHEX_TILESIZE 18
1934 /* Vector for side of triangle - ratio is close to sqrt(3) */
1935 #define GREATHEX_A 15
1936 #define GREATHEX_B 26
1937
1938 static void grid_size_greathexagonal(int width, int height,
1939                           int *tilesize, int *xextent, int *yextent)
1940 {
1941     int a = GREATHEX_A;
1942     int b = GREATHEX_B;
1943
1944     *tilesize = GREATHEX_TILESIZE;
1945     *xextent = (3*a + b) * (width-1) + 4*a;
1946     *yextent = (2*a + 2*b) * (height-1) + 3*b + a;
1947 }
1948
1949 static grid *grid_new_greathexagonal(int width, int height, const char *desc)
1950 {
1951     int x, y;
1952     int a = GREATHEX_A;
1953     int b = GREATHEX_B;
1954
1955     /* Upper bounds - don't have to be exact */
1956     int max_faces = 6 * (width + 1) * (height + 1);
1957     int max_dots = 6 * width * height;
1958
1959     tree234 *points;
1960
1961     grid *g = grid_empty();
1962     g->tilesize = GREATHEX_TILESIZE;
1963     g->faces = snewn(max_faces, grid_face);
1964     g->dots = snewn(max_dots, grid_dot);
1965
1966     points = newtree234(grid_point_cmp_fn);
1967
1968     for (y = 0; y < height; y++) {
1969         for (x = 0; x < width; x++) {
1970             grid_dot *d;
1971             /* centre of hexagon */
1972             int px = (3*a + b) * x;
1973             int py = (2*a + 2*b) * y;
1974             if (x % 2)
1975                 py += a + b;
1976
1977             /* hexagon */
1978             grid_face_add_new(g, 6);
1979             d = grid_get_dot(g, points, px - a, py - b);
1980             grid_face_set_dot(g, d, 0);
1981             d = grid_get_dot(g, points, px + a, py - b);
1982             grid_face_set_dot(g, d, 1);
1983             d = grid_get_dot(g, points, px + 2*a, py);
1984             grid_face_set_dot(g, d, 2);
1985             d = grid_get_dot(g, points, px + a, py + b);
1986             grid_face_set_dot(g, d, 3);
1987             d = grid_get_dot(g, points, px - a, py + b);
1988             grid_face_set_dot(g, d, 4);
1989             d = grid_get_dot(g, points, px - 2*a, py);
1990             grid_face_set_dot(g, d, 5);
1991
1992             /* square below hexagon */
1993             if (y < height - 1) {
1994                 grid_face_add_new(g, 4);
1995                 d = grid_get_dot(g, points, px - a, py + b);
1996                 grid_face_set_dot(g, d, 0);
1997                 d = grid_get_dot(g, points, px + a, py + b);
1998                 grid_face_set_dot(g, d, 1);
1999                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
2000                 grid_face_set_dot(g, d, 2);
2001                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
2002                 grid_face_set_dot(g, d, 3);
2003             }
2004
2005             /* square below right */
2006             if ((x < width - 1) && (((x % 2) == 0) || (y < height - 1))) {
2007                 grid_face_add_new(g, 4);
2008                 d = grid_get_dot(g, points, px + 2*a, py);
2009                 grid_face_set_dot(g, d, 0);
2010                 d = grid_get_dot(g, points, px + 2*a + b, py + a);
2011                 grid_face_set_dot(g, d, 1);
2012                 d = grid_get_dot(g, points, px + a + b, py + a + b);
2013                 grid_face_set_dot(g, d, 2);
2014                 d = grid_get_dot(g, points, px + a, py + b);
2015                 grid_face_set_dot(g, d, 3);
2016             }
2017
2018             /* square below left */
2019             if ((x > 0) && (((x % 2) == 0) || (y < height - 1))) {
2020                 grid_face_add_new(g, 4);
2021                 d = grid_get_dot(g, points, px - 2*a, py);
2022                 grid_face_set_dot(g, d, 0);
2023                 d = grid_get_dot(g, points, px - a, py + b);
2024                 grid_face_set_dot(g, d, 1);
2025                 d = grid_get_dot(g, points, px - a - b, py + a + b);
2026                 grid_face_set_dot(g, d, 2);
2027                 d = grid_get_dot(g, points, px - 2*a - b, py + a);
2028                 grid_face_set_dot(g, d, 3);
2029             }
2030
2031             /* Triangle below right */
2032             if ((x < width - 1) && (y < height - 1)) {
2033                 grid_face_add_new(g, 3);
2034                 d = grid_get_dot(g, points, px + a, py + b);
2035                 grid_face_set_dot(g, d, 0);
2036                 d = grid_get_dot(g, points, px + a + b, py + a + b);
2037                 grid_face_set_dot(g, d, 1);
2038                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
2039                 grid_face_set_dot(g, d, 2);
2040             }
2041
2042             /* Triangle below left */
2043             if ((x > 0) && (y < height - 1)) {
2044                 grid_face_add_new(g, 3);
2045                 d = grid_get_dot(g, points, px - a, py + b);
2046                 grid_face_set_dot(g, d, 0);
2047                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
2048                 grid_face_set_dot(g, d, 1);
2049                 d = grid_get_dot(g, points, px - a - b, py + a + b);
2050                 grid_face_set_dot(g, d, 2);
2051             }
2052         }
2053     }
2054
2055     freetree234(points);
2056     assert(g->num_faces <= max_faces);
2057     assert(g->num_dots <= max_dots);
2058
2059     grid_make_consistent(g);
2060     return g;
2061 }
2062
2063 #define OCTAGONAL_TILESIZE 40
2064 /* b/a approx sqrt(2) */
2065 #define OCTAGONAL_A 29
2066 #define OCTAGONAL_B 41
2067
2068 static void grid_size_octagonal(int width, int height,
2069                           int *tilesize, int *xextent, int *yextent)
2070 {
2071     int a = OCTAGONAL_A;
2072     int b = OCTAGONAL_B;
2073
2074     *tilesize = OCTAGONAL_TILESIZE;
2075     *xextent = (2*a + b) * width;
2076     *yextent = (2*a + b) * height;
2077 }
2078
2079 static grid *grid_new_octagonal(int width, int height, const char *desc)
2080 {
2081     int x, y;
2082     int a = OCTAGONAL_A;
2083     int b = OCTAGONAL_B;
2084
2085     /* Upper bounds - don't have to be exact */
2086     int max_faces = 2 * width * height;
2087     int max_dots = 4 * (width + 1) * (height + 1);
2088
2089     tree234 *points;
2090
2091     grid *g = grid_empty();
2092     g->tilesize = OCTAGONAL_TILESIZE;
2093     g->faces = snewn(max_faces, grid_face);
2094     g->dots = snewn(max_dots, grid_dot);
2095
2096     points = newtree234(grid_point_cmp_fn);
2097
2098     for (y = 0; y < height; y++) {
2099         for (x = 0; x < width; x++) {
2100             grid_dot *d;
2101             /* cell position */
2102             int px = (2*a + b) * x;
2103             int py = (2*a + b) * y;
2104             /* octagon */
2105             grid_face_add_new(g, 8);
2106             d = grid_get_dot(g, points, px + a, py);
2107             grid_face_set_dot(g, d, 0);
2108             d = grid_get_dot(g, points, px + a + b, py);
2109             grid_face_set_dot(g, d, 1);
2110             d = grid_get_dot(g, points, px + 2*a + b, py + a);
2111             grid_face_set_dot(g, d, 2);
2112             d = grid_get_dot(g, points, px + 2*a + b, py + a + b);
2113             grid_face_set_dot(g, d, 3);
2114             d = grid_get_dot(g, points, px + a + b, py + 2*a + b);
2115             grid_face_set_dot(g, d, 4);
2116             d = grid_get_dot(g, points, px + a, py + 2*a + b);
2117             grid_face_set_dot(g, d, 5);
2118             d = grid_get_dot(g, points, px, py + a + b);
2119             grid_face_set_dot(g, d, 6);
2120             d = grid_get_dot(g, points, px, py + a);
2121             grid_face_set_dot(g, d, 7);
2122
2123             /* diamond */
2124             if ((x > 0) && (y > 0)) {
2125                 grid_face_add_new(g, 4);
2126                 d = grid_get_dot(g, points, px, py - a);
2127                 grid_face_set_dot(g, d, 0);
2128                 d = grid_get_dot(g, points, px + a, py);
2129                 grid_face_set_dot(g, d, 1);
2130                 d = grid_get_dot(g, points, px, py + a);
2131                 grid_face_set_dot(g, d, 2);
2132                 d = grid_get_dot(g, points, px - a, py);
2133                 grid_face_set_dot(g, d, 3);
2134             }
2135         }
2136     }
2137
2138     freetree234(points);
2139     assert(g->num_faces <= max_faces);
2140     assert(g->num_dots <= max_dots);
2141
2142     grid_make_consistent(g);
2143     return g;
2144 }
2145
2146 #define KITE_TILESIZE 40
2147 /* b/a approx sqrt(3) */
2148 #define KITE_A 15
2149 #define KITE_B 26
2150
2151 static void grid_size_kites(int width, int height,
2152                      int *tilesize, int *xextent, int *yextent)
2153 {
2154     int a = KITE_A;
2155     int b = KITE_B;
2156
2157     *tilesize = KITE_TILESIZE;
2158     *xextent = 4*b * width + 2*b;
2159     *yextent = 6*a * (height-1) + 8*a;
2160 }
2161
2162 static grid *grid_new_kites(int width, int height, const char *desc)
2163 {
2164     int x, y;
2165     int a = KITE_A;
2166     int b = KITE_B;
2167
2168     /* Upper bounds - don't have to be exact */
2169     int max_faces = 6 * width * height;
2170     int max_dots = 6 * (width + 1) * (height + 1);
2171
2172     tree234 *points;
2173
2174     grid *g = grid_empty();
2175     g->tilesize = KITE_TILESIZE;
2176     g->faces = snewn(max_faces, grid_face);
2177     g->dots = snewn(max_dots, grid_dot);
2178
2179     points = newtree234(grid_point_cmp_fn);
2180
2181     for (y = 0; y < height; y++) {
2182         for (x = 0; x < width; x++) {
2183             grid_dot *d;
2184             /* position of order-6 dot */
2185             int px = 4*b * x;
2186             int py = 6*a * y;
2187             if (y % 2)
2188                 px += 2*b;
2189
2190             /* kite pointing up-left */
2191             grid_face_add_new(g, 4);
2192             d = grid_get_dot(g, points, px, py);
2193             grid_face_set_dot(g, d, 0);
2194             d = grid_get_dot(g, points, px + 2*b, py);
2195             grid_face_set_dot(g, d, 1);
2196             d = grid_get_dot(g, points, px + 2*b, py + 2*a);
2197             grid_face_set_dot(g, d, 2);
2198             d = grid_get_dot(g, points, px + b, py + 3*a);
2199             grid_face_set_dot(g, d, 3);
2200
2201             /* kite pointing up */
2202             grid_face_add_new(g, 4);
2203             d = grid_get_dot(g, points, px, py);
2204             grid_face_set_dot(g, d, 0);
2205             d = grid_get_dot(g, points, px + b, py + 3*a);
2206             grid_face_set_dot(g, d, 1);
2207             d = grid_get_dot(g, points, px, py + 4*a);
2208             grid_face_set_dot(g, d, 2);
2209             d = grid_get_dot(g, points, px - b, py + 3*a);
2210             grid_face_set_dot(g, d, 3);
2211
2212             /* kite pointing up-right */
2213             grid_face_add_new(g, 4);
2214             d = grid_get_dot(g, points, px, py);
2215             grid_face_set_dot(g, d, 0);
2216             d = grid_get_dot(g, points, px - b, py + 3*a);
2217             grid_face_set_dot(g, d, 1);
2218             d = grid_get_dot(g, points, px - 2*b, py + 2*a);
2219             grid_face_set_dot(g, d, 2);
2220             d = grid_get_dot(g, points, px - 2*b, py);
2221             grid_face_set_dot(g, d, 3);
2222
2223             /* kite pointing down-right */
2224             grid_face_add_new(g, 4);
2225             d = grid_get_dot(g, points, px, py);
2226             grid_face_set_dot(g, d, 0);
2227             d = grid_get_dot(g, points, px - 2*b, py);
2228             grid_face_set_dot(g, d, 1);
2229             d = grid_get_dot(g, points, px - 2*b, py - 2*a);
2230             grid_face_set_dot(g, d, 2);
2231             d = grid_get_dot(g, points, px - b, py - 3*a);
2232             grid_face_set_dot(g, d, 3);
2233
2234             /* kite pointing down */
2235             grid_face_add_new(g, 4);
2236             d = grid_get_dot(g, points, px, py);
2237             grid_face_set_dot(g, d, 0);
2238             d = grid_get_dot(g, points, px - b, py - 3*a);
2239             grid_face_set_dot(g, d, 1);
2240             d = grid_get_dot(g, points, px, py - 4*a);
2241             grid_face_set_dot(g, d, 2);
2242             d = grid_get_dot(g, points, px + b, py - 3*a);
2243             grid_face_set_dot(g, d, 3);
2244
2245             /* kite pointing down-left */
2246             grid_face_add_new(g, 4);
2247             d = grid_get_dot(g, points, px, py);
2248             grid_face_set_dot(g, d, 0);
2249             d = grid_get_dot(g, points, px + b, py - 3*a);
2250             grid_face_set_dot(g, d, 1);
2251             d = grid_get_dot(g, points, px + 2*b, py - 2*a);
2252             grid_face_set_dot(g, d, 2);
2253             d = grid_get_dot(g, points, px + 2*b, py);
2254             grid_face_set_dot(g, d, 3);
2255         }
2256     }
2257
2258     freetree234(points);
2259     assert(g->num_faces <= max_faces);
2260     assert(g->num_dots <= max_dots);
2261
2262     grid_make_consistent(g);
2263     return g;
2264 }
2265
2266 #define FLORET_TILESIZE 150
2267 /* -py/px is close to tan(30 - atan(sqrt(3)/9))
2268  * using py=26 makes everything lean to the left, rather than right
2269  */
2270 #define FLORET_PX 75
2271 #define FLORET_PY -26
2272
2273 static void grid_size_floret(int width, int height,
2274                           int *tilesize, int *xextent, int *yextent)
2275 {
2276     int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
2277     int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
2278     int ry = qy-py;
2279     /* rx unused in determining grid size. */
2280
2281     *tilesize = FLORET_TILESIZE;
2282     *xextent = (6*px+3*qx)/2 * (width-1) + 4*qx + 2*px;
2283     *yextent = (5*qy-4*py) * (height-1) + 4*qy + 2*ry;
2284 }
2285
2286 static grid *grid_new_floret(int width, int height, const char *desc)
2287 {
2288     int x, y;
2289     /* Vectors for sides; weird numbers needed to keep puzzle aligned with window
2290      * -py/px is close to tan(30 - atan(sqrt(3)/9))
2291      * using py=26 makes everything lean to the left, rather than right
2292      */
2293     int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
2294     int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
2295     int rx = qx-px, ry = qy-py;                 /* |(-15,  78)| = 79.38 */
2296
2297     /* Upper bounds - don't have to be exact */
2298     int max_faces = 6 * width * height;
2299     int max_dots = 9 * (width + 1) * (height + 1);
2300
2301     tree234 *points;
2302
2303     grid *g = grid_empty();
2304     g->tilesize = FLORET_TILESIZE;
2305     g->faces = snewn(max_faces, grid_face);
2306     g->dots = snewn(max_dots, grid_dot);
2307
2308     points = newtree234(grid_point_cmp_fn);
2309
2310     /* generate pentagonal faces */
2311     for (y = 0; y < height; y++) {
2312         for (x = 0; x < width; x++) {
2313             grid_dot *d;
2314             /* face centre */
2315             int cx = (6*px+3*qx)/2 * x;
2316             int cy = (4*py-5*qy) * y;
2317             if (x % 2)
2318                 cy -= (4*py-5*qy)/2;
2319             else if (y && y == height-1)
2320                 continue; /* make better looking grids?  try 3x3 for instance */
2321
2322             grid_face_add_new(g, 5);
2323             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2324             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 1);
2325             d = grid_get_dot(g, points, cx+2*rx+qx, cy+2*ry+qy); grid_face_set_dot(g, d, 2);
2326             d = grid_get_dot(g, points, cx+2*qx+rx, cy+2*qy+ry); grid_face_set_dot(g, d, 3);
2327             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 4);
2328
2329             grid_face_add_new(g, 5);
2330             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2331             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 1);
2332             d = grid_get_dot(g, points, cx+2*qx+px, cy+2*qy+py); grid_face_set_dot(g, d, 2);
2333             d = grid_get_dot(g, points, cx+2*px+qx, cy+2*py+qy); grid_face_set_dot(g, d, 3);
2334             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 4);
2335
2336             grid_face_add_new(g, 5);
2337             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2338             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 1);
2339             d = grid_get_dot(g, points, cx+2*px-rx, cy+2*py-ry); grid_face_set_dot(g, d, 2);
2340             d = grid_get_dot(g, points, cx-2*rx+px, cy-2*ry+py); grid_face_set_dot(g, d, 3);
2341             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 4);
2342
2343             grid_face_add_new(g, 5);
2344             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2345             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 1);
2346             d = grid_get_dot(g, points, cx-2*rx-qx, cy-2*ry-qy); grid_face_set_dot(g, d, 2);
2347             d = grid_get_dot(g, points, cx-2*qx-rx, cy-2*qy-ry); grid_face_set_dot(g, d, 3);
2348             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 4);
2349
2350             grid_face_add_new(g, 5);
2351             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2352             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 1);
2353             d = grid_get_dot(g, points, cx-2*qx-px, cy-2*qy-py); grid_face_set_dot(g, d, 2);
2354             d = grid_get_dot(g, points, cx-2*px-qx, cy-2*py-qy); grid_face_set_dot(g, d, 3);
2355             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 4);
2356
2357             grid_face_add_new(g, 5);
2358             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2359             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 1);
2360             d = grid_get_dot(g, points, cx-2*px+rx, cy-2*py+ry); grid_face_set_dot(g, d, 2);
2361             d = grid_get_dot(g, points, cx+2*rx-px, cy+2*ry-py); grid_face_set_dot(g, d, 3);
2362             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 4);
2363         }
2364     }
2365
2366     freetree234(points);
2367     assert(g->num_faces <= max_faces);
2368     assert(g->num_dots <= max_dots);
2369
2370     grid_make_consistent(g);
2371     return g;
2372 }
2373
2374 /* DODEC_* are used for dodecagonal and great-dodecagonal grids. */
2375 #define DODEC_TILESIZE 26
2376 /* Vector for side of triangle - ratio is close to sqrt(3) */
2377 #define DODEC_A 15
2378 #define DODEC_B 26
2379
2380 static void grid_size_dodecagonal(int width, int height,
2381                           int *tilesize, int *xextent, int *yextent)
2382 {
2383     int a = DODEC_A;
2384     int b = DODEC_B;
2385
2386     *tilesize = DODEC_TILESIZE;
2387     *xextent = (4*a + 2*b) * (width-1) + 3*(2*a + b);
2388     *yextent = (3*a + 2*b) * (height-1) + 2*(2*a + b);
2389 }
2390
2391 static grid *grid_new_dodecagonal(int width, int height, const char *desc)
2392 {
2393     int x, y;
2394     int a = DODEC_A;
2395     int b = DODEC_B;
2396
2397     /* Upper bounds - don't have to be exact */
2398     int max_faces = 3 * width * height;
2399     int max_dots = 14 * width * height;
2400
2401     tree234 *points;
2402
2403     grid *g = grid_empty();
2404     g->tilesize = DODEC_TILESIZE;
2405     g->faces = snewn(max_faces, grid_face);
2406     g->dots = snewn(max_dots, grid_dot);
2407
2408     points = newtree234(grid_point_cmp_fn);
2409
2410     for (y = 0; y < height; y++) {
2411         for (x = 0; x < width; x++) {
2412             grid_dot *d;
2413             /* centre of dodecagon */
2414             int px = (4*a + 2*b) * x;
2415             int py = (3*a + 2*b) * y;
2416             if (y % 2)
2417                 px += 2*a + b;
2418
2419             /* dodecagon */
2420             grid_face_add_new(g, 12);
2421             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
2422             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
2423             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
2424             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
2425             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
2426             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
2427             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
2428             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
2429             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
2430             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
2431             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
2432             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
2433
2434             /* triangle below dodecagon */
2435             if ((y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
2436                 grid_face_add_new(g, 3);
2437                 d = grid_get_dot(g, points, px + a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2438                 d = grid_get_dot(g, points, px    , py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2439                 d = grid_get_dot(g, points, px - a, py + (2*a +   b)); grid_face_set_dot(g, d, 2);
2440             }
2441
2442             /* triangle above dodecagon */
2443             if ((y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
2444                 grid_face_add_new(g, 3);
2445                 d = grid_get_dot(g, points, px - a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2446                 d = grid_get_dot(g, points, px    , py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2447                 d = grid_get_dot(g, points, px + a, py - (2*a +   b)); grid_face_set_dot(g, d, 2);
2448             }
2449         }
2450     }
2451
2452     freetree234(points);
2453     assert(g->num_faces <= max_faces);
2454     assert(g->num_dots <= max_dots);
2455
2456     grid_make_consistent(g);
2457     return g;
2458 }
2459
2460 static void grid_size_greatdodecagonal(int width, int height,
2461                           int *tilesize, int *xextent, int *yextent)
2462 {
2463     int a = DODEC_A;
2464     int b = DODEC_B;
2465
2466     *tilesize = DODEC_TILESIZE;
2467     *xextent = (6*a + 2*b) * (width-1) + 2*(2*a + b) + 3*a + b;
2468     *yextent = (3*a + 3*b) * (height-1) + 2*(2*a + b);
2469 }
2470
2471 static grid *grid_new_greatdodecagonal(int width, int height, const char *desc)
2472 {
2473     int x, y;
2474     /* Vector for side of triangle - ratio is close to sqrt(3) */
2475     int a = DODEC_A;
2476     int b = DODEC_B;
2477
2478     /* Upper bounds - don't have to be exact */
2479     int max_faces = 30 * width * height;
2480     int max_dots = 200 * width * height;
2481
2482     tree234 *points;
2483
2484     grid *g = grid_empty();
2485     g->tilesize = DODEC_TILESIZE;
2486     g->faces = snewn(max_faces, grid_face);
2487     g->dots = snewn(max_dots, grid_dot);
2488
2489     points = newtree234(grid_point_cmp_fn);
2490
2491     for (y = 0; y < height; y++) {
2492         for (x = 0; x < width; x++) {
2493             grid_dot *d;
2494             /* centre of dodecagon */
2495             int px = (6*a + 2*b) * x;
2496             int py = (3*a + 3*b) * y;
2497             if (y % 2)
2498                 px += 3*a + b;
2499
2500             /* dodecagon */
2501             grid_face_add_new(g, 12);
2502             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
2503             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
2504             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
2505             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
2506             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
2507             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
2508             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
2509             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
2510             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
2511             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
2512             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
2513             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
2514
2515             /* hexagon below dodecagon */
2516             if (y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2517                 grid_face_add_new(g, 6);
2518                 d = grid_get_dot(g, points, px +   a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2519                 d = grid_get_dot(g, points, px + 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2520                 d = grid_get_dot(g, points, px +   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2521                 d = grid_get_dot(g, points, px -   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2522                 d = grid_get_dot(g, points, px - 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2523                 d = grid_get_dot(g, points, px -   a, py + (2*a +   b)); grid_face_set_dot(g, d, 5);
2524             }
2525
2526             /* hexagon above dodecagon */
2527             if (y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2528                 grid_face_add_new(g, 6);
2529                 d = grid_get_dot(g, points, px -   a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2530                 d = grid_get_dot(g, points, px - 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2531                 d = grid_get_dot(g, points, px -   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2532                 d = grid_get_dot(g, points, px +   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2533                 d = grid_get_dot(g, points, px + 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2534                 d = grid_get_dot(g, points, px +   a, py - (2*a +   b)); grid_face_set_dot(g, d, 5);
2535             }
2536
2537             /* square on right of dodecagon */
2538             if (x < width - 1) {
2539                 grid_face_add_new(g, 4);
2540                 d = grid_get_dot(g, points, px + 2*a + b, py - a); grid_face_set_dot(g, d, 0);
2541                 d = grid_get_dot(g, points, px + 4*a + b, py - a); grid_face_set_dot(g, d, 1);
2542                 d = grid_get_dot(g, points, px + 4*a + b, py + a); grid_face_set_dot(g, d, 2);
2543                 d = grid_get_dot(g, points, px + 2*a + b, py + a); grid_face_set_dot(g, d, 3);
2544             }
2545
2546             /* square on top right of dodecagon */
2547             if (y && (x < width - 1 || !(y % 2))) {
2548                 grid_face_add_new(g, 4);
2549                 d = grid_get_dot(g, points, px + (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2550                 d = grid_get_dot(g, points, px + (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2551                 d = grid_get_dot(g, points, px + (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 2);
2552                 d = grid_get_dot(g, points, px + (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 3);
2553             }
2554
2555             /* square on top left of dodecagon */
2556             if (y && (x || (y % 2))) {
2557                 grid_face_add_new(g, 4);
2558                 d = grid_get_dot(g, points, px - (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 0);
2559                 d = grid_get_dot(g, points, px - (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 1);
2560                 d = grid_get_dot(g, points, px - (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 2);
2561                 d = grid_get_dot(g, points, px - (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 3);
2562             }
2563         }
2564     }
2565
2566     freetree234(points);
2567     assert(g->num_faces <= max_faces);
2568     assert(g->num_dots <= max_dots);
2569
2570     grid_make_consistent(g);
2571     return g;
2572 }
2573
2574 static void grid_size_greatgreatdodecagonal(int width, int height,
2575                                             int *tilesize, int *xextent, int *yextent)
2576 {
2577     int a = DODEC_A;
2578     int b = DODEC_B;
2579
2580     *tilesize = DODEC_TILESIZE;
2581     *xextent = (4*a + 4*b) * (width-1) + 2*(2*a + b) + 2*a + 2*b;
2582     *yextent = (6*a + 2*b) * (height-1) + 2*(2*a + b);
2583 }
2584
2585 static grid *grid_new_greatgreatdodecagonal(int width, int height, const char *desc)
2586 {
2587     int x, y;
2588     /* Vector for side of triangle - ratio is close to sqrt(3) */
2589     int a = DODEC_A;
2590     int b = DODEC_B;
2591
2592     /* Upper bounds - don't have to be exact */
2593     int max_faces = 50 * width * height;
2594     int max_dots = 300 * width * height;
2595
2596     tree234 *points;
2597
2598     grid *g = grid_empty();
2599     g->tilesize = DODEC_TILESIZE;
2600     g->faces = snewn(max_faces, grid_face);
2601     g->dots = snewn(max_dots, grid_dot);
2602
2603     points = newtree234(grid_point_cmp_fn);
2604
2605     for (y = 0; y < height; y++) {
2606         for (x = 0; x < width; x++) {
2607             grid_dot *d;
2608             /* centre of dodecagon */
2609             int px = (4*a + 4*b) * x;
2610             int py = (6*a + 2*b) * y;
2611             if (y % 2)
2612                 px += 2*a + 2*b;
2613
2614             /* dodecagon */
2615             grid_face_add_new(g, 12);
2616             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
2617             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
2618             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
2619             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
2620             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
2621             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
2622             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
2623             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
2624             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
2625             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
2626             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
2627             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
2628
2629             /* hexagon on top right of dodecagon */
2630             if (y && (x < width - 1 || !(y % 2))) {
2631                 grid_face_add_new(g, 6);
2632                 d = grid_get_dot(g, points, px + (a + 2*b), py - (4*a + b)); grid_face_set_dot(g, d, 0);
2633                 d = grid_get_dot(g, points, px + (a + 2*b), py - (2*a + b)); grid_face_set_dot(g, d, 1);
2634                 d = grid_get_dot(g, points, px + (a +   b), py - (  a + b)); grid_face_set_dot(g, d, 2);
2635                 d = grid_get_dot(g, points, px + (a      ), py - (2*a + b)); grid_face_set_dot(g, d, 3);
2636                 d = grid_get_dot(g, points, px + (a      ), py - (4*a + b)); grid_face_set_dot(g, d, 4);
2637                 d = grid_get_dot(g, points, px + (a +   b), py - (5*a + b)); grid_face_set_dot(g, d, 5);
2638             }
2639
2640             /* hexagon on right of dodecagon*/
2641             if (x < width - 1) {
2642                 grid_face_add_new(g, 6);
2643                 d = grid_get_dot(g, points, px + (2*a + 3*b), py -   a); grid_face_set_dot(g, d, 0);
2644                 d = grid_get_dot(g, points, px + (2*a + 3*b), py +   a); grid_face_set_dot(g, d, 1);
2645                 d = grid_get_dot(g, points, px + (2*a + 2*b), py + 2*a); grid_face_set_dot(g, d, 2);
2646                 d = grid_get_dot(g, points, px + (2*a +   b), py +   a); grid_face_set_dot(g, d, 3);
2647                 d = grid_get_dot(g, points, px + (2*a +   b), py -   a); grid_face_set_dot(g, d, 4);
2648                 d = grid_get_dot(g, points, px + (2*a + 2*b), py - 2*a); grid_face_set_dot(g, d, 5);
2649             }
2650
2651             /* hexagon on bottom right  of dodecagon */
2652             if ((y < height - 1) && (x < width - 1 || !(y % 2))) {
2653                 grid_face_add_new(g, 6);
2654                 d = grid_get_dot(g, points, px + (a + 2*b), py + (2*a + b)); grid_face_set_dot(g, d, 0);
2655                 d = grid_get_dot(g, points, px + (a + 2*b), py + (4*a + b)); grid_face_set_dot(g, d, 1);
2656                 d = grid_get_dot(g, points, px + (a +   b), py + (5*a + b)); grid_face_set_dot(g, d, 2);
2657                 d = grid_get_dot(g, points, px + (a      ), py + (4*a + b)); grid_face_set_dot(g, d, 3);
2658                 d = grid_get_dot(g, points, px + (a      ), py + (2*a + b)); grid_face_set_dot(g, d, 4);
2659                 d = grid_get_dot(g, points, px + (a +   b), py + (  a + b)); grid_face_set_dot(g, d, 5);
2660             }
2661
2662             /* square on top right of dodecagon */
2663             if (y && (x < width - 1 )) {
2664                 grid_face_add_new(g, 4);
2665                 d = grid_get_dot(g, points, px + (  a + 2*b), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2666                 d = grid_get_dot(g, points, px + (2*a + 2*b), py - (2*a      )); grid_face_set_dot(g, d, 1);
2667                 d = grid_get_dot(g, points, px + (2*a +   b), py - (  a      )); grid_face_set_dot(g, d, 2);
2668                 d = grid_get_dot(g, points, px + (  a +   b), py - (  a +   b)); grid_face_set_dot(g, d, 3);
2669             }
2670
2671             /* square on bottom right of dodecagon */
2672             if ((y < height - 1) && (x < width - 1 )) {
2673                 grid_face_add_new(g, 4);
2674                 d = grid_get_dot(g, points, px + (2*a + 2*b), py + (2*a      )); grid_face_set_dot(g, d, 0);
2675                 d = grid_get_dot(g, points, px + (  a + 2*b), py + (2*a +   b)); grid_face_set_dot(g, d, 1);
2676                 d = grid_get_dot(g, points, px + (  a +   b), py + (  a +   b)); grid_face_set_dot(g, d, 2);
2677                 d = grid_get_dot(g, points, px + (2*a +   b), py + (  a      )); grid_face_set_dot(g, d, 3);
2678             }
2679
2680             /* square below dodecagon */
2681             if ((y < height - 1) && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2682                 grid_face_add_new(g, 4);
2683                 d = grid_get_dot(g, points, px + a, py + (2*a + b)); grid_face_set_dot(g, d, 0);
2684                 d = grid_get_dot(g, points, px + a, py + (4*a + b)); grid_face_set_dot(g, d, 1);
2685                 d = grid_get_dot(g, points, px - a, py + (4*a + b)); grid_face_set_dot(g, d, 2);
2686                 d = grid_get_dot(g, points, px - a, py + (2*a + b)); grid_face_set_dot(g, d, 3);
2687             }
2688
2689             /* square on bottom left of dodecagon */
2690             if (x && (y < height - 1)) {
2691                 grid_face_add_new(g, 4);
2692                 d = grid_get_dot(g, points, px - (2*a +   b), py + (  a      )); grid_face_set_dot(g, d, 0);
2693                 d = grid_get_dot(g, points, px - (  a +   b), py + (  a +   b)); grid_face_set_dot(g, d, 1);
2694                 d = grid_get_dot(g, points, px - (  a + 2*b), py + (2*a +   b)); grid_face_set_dot(g, d, 2);
2695                 d = grid_get_dot(g, points, px - (2*a + 2*b), py + (2*a      )); grid_face_set_dot(g, d, 3);
2696             }
2697
2698             /* square on top left of dodecagon */
2699             if (x && y) {
2700                 grid_face_add_new(g, 4);
2701                 d = grid_get_dot(g, points, px - (  a +   b), py - (  a +   b)); grid_face_set_dot(g, d, 0);
2702                 d = grid_get_dot(g, points, px - (2*a +   b), py - (  a      )); grid_face_set_dot(g, d, 1);
2703                 d = grid_get_dot(g, points, px - (2*a + 2*b), py - (2*a      )); grid_face_set_dot(g, d, 2);
2704                 d = grid_get_dot(g, points, px - (  a + 2*b), py - (2*a +   b)); grid_face_set_dot(g, d, 3);
2705
2706             }
2707
2708             /* square above dodecagon */
2709             if (y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2710                 grid_face_add_new(g, 4);
2711                 d = grid_get_dot(g, points, px + a, py - (4*a + b)); grid_face_set_dot(g, d, 0);
2712                 d = grid_get_dot(g, points, px + a, py - (2*a + b)); grid_face_set_dot(g, d, 1);
2713                 d = grid_get_dot(g, points, px - a, py - (2*a + b)); grid_face_set_dot(g, d, 2);
2714                 d = grid_get_dot(g, points, px - a, py - (4*a + b)); grid_face_set_dot(g, d, 3);
2715             }
2716
2717             /* upper triangle (v) */
2718             if (y && (x < width - 1)) {
2719                 grid_face_add_new(g, 3);
2720                 d = grid_get_dot(g, points, px + (3*a + 2*b), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2721                 d = grid_get_dot(g, points, px + (2*a + 2*b), py - (2*a      )); grid_face_set_dot(g, d, 1);
2722                 d = grid_get_dot(g, points, px + (  a + 2*b), py - (2*a +   b)); grid_face_set_dot(g, d, 2);
2723             }
2724
2725             /* lower triangle (^) */
2726             if ((y < height - 1) && (x < width - 1)) {
2727                 grid_face_add_new(g, 3);
2728                 d = grid_get_dot(g, points, px + (3*a + 2*b), py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2729                 d = grid_get_dot(g, points, px + (  a + 2*b), py + (2*a +   b)); grid_face_set_dot(g, d, 1);
2730                 d = grid_get_dot(g, points, px + (2*a + 2*b), py + (2*a      )); grid_face_set_dot(g, d, 2);
2731             }
2732         }
2733     }
2734
2735     freetree234(points);
2736     assert(g->num_faces <= max_faces);
2737     assert(g->num_dots <= max_dots);
2738
2739     grid_make_consistent(g);
2740     return g;
2741 }
2742
2743 typedef struct setface_ctx
2744 {
2745     int xmin, xmax, ymin, ymax;
2746
2747     grid *g;
2748     tree234 *points;
2749 } setface_ctx;
2750
2751 static double round_int_nearest_away(double r)
2752 {
2753     return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
2754 }
2755
2756 static int set_faces(penrose_state *state, vector *vs, int n, int depth)
2757 {
2758     setface_ctx *sf_ctx = (setface_ctx *)state->ctx;
2759     int i;
2760     int xs[4], ys[4];
2761
2762     if (depth < state->max_depth) return 0;
2763 #ifdef DEBUG_PENROSE
2764     if (n != 4) return 0; /* triangles are sent as debugging. */
2765 #endif
2766
2767     for (i = 0; i < n; i++) {
2768         double tx = v_x(vs, i), ty = v_y(vs, i);
2769
2770         xs[i] = (int)round_int_nearest_away(tx);
2771         ys[i] = (int)round_int_nearest_away(ty);
2772
2773         if (xs[i] < sf_ctx->xmin || xs[i] > sf_ctx->xmax) return 0;
2774         if (ys[i] < sf_ctx->ymin || ys[i] > sf_ctx->ymax) return 0;
2775     }
2776
2777     grid_face_add_new(sf_ctx->g, n);
2778     debug(("penrose: new face l=%f gen=%d...",
2779            penrose_side_length(state->start_size, depth), depth));
2780     for (i = 0; i < n; i++) {
2781         grid_dot *d = grid_get_dot(sf_ctx->g, sf_ctx->points,
2782                                    xs[i], ys[i]);
2783         grid_face_set_dot(sf_ctx->g, d, i);
2784         debug((" ... dot 0x%x (%d,%d) (was %2.2f,%2.2f)",
2785                d, d->x, d->y, v_x(vs, i), v_y(vs, i)));
2786     }
2787
2788     return 0;
2789 }
2790
2791 #define PENROSE_TILESIZE 100
2792
2793 static void grid_size_penrose(int width, int height,
2794                        int *tilesize, int *xextent, int *yextent)
2795 {
2796     int l = PENROSE_TILESIZE;
2797
2798     *tilesize = l;
2799     *xextent = l * width;
2800     *yextent = l * height;
2801 }
2802
2803 static grid *grid_new_penrose(int width, int height, int which, const char *desc); /* forward reference */
2804
2805 static char *grid_new_desc_penrose(grid_type type, int width, int height, random_state *rs)
2806 {
2807     int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff;
2808     double outer_radius;
2809     int inner_radius;
2810     char gd[255];
2811     int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
2812     grid *g;
2813
2814     while (1) {
2815         /* We want to produce a random bit of penrose tiling, so we
2816          * calculate a random offset (within the patch that penrose.c
2817          * calculates for us) and an angle (multiple of 36) to rotate
2818          * the patch. */
2819
2820         penrose_calculate_size(which, tilesize, width, height,
2821                                &outer_radius, &startsz, &depth);
2822
2823         /* Calculate radius of (circumcircle of) patch, subtract from
2824          * radius calculated. */
2825         inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
2826
2827         /* Pick a random offset (the easy way: choose within outer
2828          * square, discarding while it's outside the circle) */
2829         do {
2830             xoff = random_upto(rs, 2*inner_radius) - inner_radius;
2831             yoff = random_upto(rs, 2*inner_radius) - inner_radius;
2832         } while (sqrt(xoff*xoff+yoff*yoff) > inner_radius);
2833
2834         aoff = random_upto(rs, 360/36) * 36;
2835
2836         debug(("grid_desc: ts %d, %dx%d patch, orad %2.2f irad %d",
2837                tilesize, width, height, outer_radius, inner_radius));
2838         debug(("    -> xoff %d yoff %d aoff %d", xoff, yoff, aoff));
2839
2840         sprintf(gd, "G%d,%d,%d", xoff, yoff, aoff);
2841
2842         /*
2843          * Now test-generate our grid, to make sure it actually
2844          * produces something.
2845          */
2846         g = grid_new_penrose(width, height, which, gd);
2847         if (g) {
2848             grid_free(g);
2849             break;
2850         }
2851         /* If not, go back to the top of this while loop and try again
2852          * with a different random offset. */
2853     }
2854
2855     return dupstr(gd);
2856 }
2857
2858 static char *grid_validate_desc_penrose(grid_type type, int width, int height,
2859                                         const char *desc)
2860 {
2861     int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff, inner_radius;
2862     double outer_radius;
2863     int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
2864     grid *g;
2865
2866     if (!desc)
2867         return "Missing grid description string.";
2868
2869     penrose_calculate_size(which, tilesize, width, height,
2870                            &outer_radius, &startsz, &depth);
2871     inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
2872
2873     if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
2874         return "Invalid format grid description string.";
2875
2876     if (sqrt(xoff*xoff + yoff*yoff) > inner_radius)
2877         return "Patch offset out of bounds.";
2878     if ((aoff % 36) != 0 || aoff < 0 || aoff >= 360)
2879         return "Angle offset out of bounds.";
2880
2881     /*
2882      * Test-generate to ensure these parameters don't end us up with
2883      * no grid at all.
2884      */
2885     g = grid_new_penrose(width, height, which, desc);
2886     if (!g)
2887         return "Patch coordinates do not identify a usable grid fragment";
2888     grid_free(g);
2889
2890     return NULL;
2891 }
2892
2893 /*
2894  * We're asked for a grid of a particular size, and we generate enough
2895  * of the tiling so we can be sure to have enough random grid from which
2896  * to pick.
2897  */
2898
2899 static grid *grid_new_penrose(int width, int height, int which, const char *desc)
2900 {
2901     int max_faces, max_dots, tilesize = PENROSE_TILESIZE;
2902     int xsz, ysz, xoff, yoff, aoff;
2903     double rradius;
2904
2905     tree234 *points;
2906     grid *g;
2907
2908     penrose_state ps;
2909     setface_ctx sf_ctx;
2910
2911     penrose_calculate_size(which, tilesize, width, height,
2912                            &rradius, &ps.start_size, &ps.max_depth);
2913
2914     debug(("penrose: w%d h%d, tile size %d, start size %d, depth %d",
2915            width, height, tilesize, ps.start_size, ps.max_depth));
2916
2917     ps.new_tile = set_faces;
2918     ps.ctx = &sf_ctx;
2919
2920     max_faces = (width*3) * (height*3); /* somewhat paranoid... */
2921     max_dots = max_faces * 4; /* ditto... */
2922
2923     g = grid_empty();
2924     g->tilesize = tilesize;
2925     g->faces = snewn(max_faces, grid_face);
2926     g->dots = snewn(max_dots, grid_dot);
2927
2928     points = newtree234(grid_point_cmp_fn);
2929
2930     memset(&sf_ctx, 0, sizeof(sf_ctx));
2931     sf_ctx.g = g;
2932     sf_ctx.points = points;
2933
2934     if (desc != NULL) {
2935         if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
2936             assert(!"Invalid grid description.");
2937     } else {
2938         xoff = yoff = aoff = 0;
2939     }
2940
2941     xsz = width * tilesize;
2942     ysz = height * tilesize;
2943
2944     sf_ctx.xmin = xoff - xsz/2;
2945     sf_ctx.xmax = xoff + xsz/2;
2946     sf_ctx.ymin = yoff - ysz/2;
2947     sf_ctx.ymax = yoff + ysz/2;
2948
2949     debug(("penrose: centre (%f, %f) xsz %f ysz %f",
2950            0.0, 0.0, xsz, ysz));
2951     debug(("penrose: x range (%f --> %f), y range (%f --> %f)",
2952            sf_ctx.xmin, sf_ctx.xmax, sf_ctx.ymin, sf_ctx.ymax));
2953
2954     penrose(&ps, which, aoff);
2955
2956     freetree234(points);
2957     assert(g->num_faces <= max_faces);
2958     assert(g->num_dots <= max_dots);
2959
2960     debug(("penrose: %d faces total (equivalent to %d wide by %d high)",
2961            g->num_faces, g->num_faces/height, g->num_faces/width));
2962
2963     /*
2964      * Return NULL if we ended up with an empty grid, either because
2965      * the initial generation was over too small a rectangle to
2966      * encompass any face or because grid_trim_vigorously ended up
2967      * removing absolutely everything.
2968      */
2969     if (g->num_faces == 0 || g->num_dots == 0) {
2970         grid_free(g);
2971         return NULL;
2972     }
2973     grid_trim_vigorously(g);
2974     if (g->num_faces == 0 || g->num_dots == 0) {
2975         grid_free(g);
2976         return NULL;
2977     }
2978
2979     grid_make_consistent(g);
2980
2981     /*
2982      * Centre the grid in its originally promised rectangle.
2983      */
2984     g->lowest_x -= ((sf_ctx.xmax - sf_ctx.xmin) -
2985                     (g->highest_x - g->lowest_x)) / 2;
2986     g->highest_x = g->lowest_x + (sf_ctx.xmax - sf_ctx.xmin);
2987     g->lowest_y -= ((sf_ctx.ymax - sf_ctx.ymin) -
2988                     (g->highest_y - g->lowest_y)) / 2;
2989     g->highest_y = g->lowest_y + (sf_ctx.ymax - sf_ctx.ymin);
2990
2991     return g;
2992 }
2993
2994 static void grid_size_penrose_p2_kite(int width, int height,
2995                        int *tilesize, int *xextent, int *yextent)
2996 {
2997     grid_size_penrose(width, height, tilesize, xextent, yextent);
2998 }
2999
3000 static void grid_size_penrose_p3_thick(int width, int height,
3001                        int *tilesize, int *xextent, int *yextent)
3002 {
3003     grid_size_penrose(width, height, tilesize, xextent, yextent);
3004 }
3005
3006 static grid *grid_new_penrose_p2_kite(int width, int height, const char *desc)
3007 {
3008     return grid_new_penrose(width, height, PENROSE_P2, desc);
3009 }
3010
3011 static grid *grid_new_penrose_p3_thick(int width, int height, const char *desc)
3012 {
3013     return grid_new_penrose(width, height, PENROSE_P3, desc);
3014 }
3015
3016 /* ----------- End of grid generators ------------- */
3017
3018 #define FNNEW(upper,lower) &grid_new_ ## lower,
3019 #define FNSZ(upper,lower) &grid_size_ ## lower,
3020
3021 static grid *(*(grid_news[]))(int, int, const char*) = { GRIDGEN_LIST(FNNEW) };
3022 static void(*(grid_sizes[]))(int, int, int*, int*, int*) = { GRIDGEN_LIST(FNSZ) };
3023
3024 char *grid_new_desc(grid_type type, int width, int height, random_state *rs)
3025 {
3026     if (type == GRID_PENROSE_P2 || type == GRID_PENROSE_P3) {
3027         return grid_new_desc_penrose(type, width, height, rs);
3028     } else if (type == GRID_TRIANGULAR) {
3029         return dupstr("0"); /* up-to-date version of triangular grid */
3030     } else {
3031         return NULL;
3032     }
3033 }
3034
3035 char *grid_validate_desc(grid_type type, int width, int height,
3036                          const char *desc)
3037 {
3038     if (type == GRID_PENROSE_P2 || type == GRID_PENROSE_P3) {
3039         return grid_validate_desc_penrose(type, width, height, desc);
3040     } else if (type == GRID_TRIANGULAR) {
3041         return grid_validate_desc_triangular(type, width, height, desc);
3042     } else {
3043         if (desc != NULL)
3044             return "Grid description strings not used with this grid type";
3045         return NULL;
3046     }
3047 }
3048
3049 grid *grid_new(grid_type type, int width, int height, const char *desc)
3050 {
3051     char *err = grid_validate_desc(type, width, height, desc);
3052     if (err) assert(!"Invalid grid description.");
3053
3054     return grid_news[type](width, height, desc);
3055 }
3056
3057 void grid_compute_size(grid_type type, int width, int height,
3058                        int *tilesize, int *xextent, int *yextent)
3059 {
3060     grid_sizes[type](width, height, tilesize, xextent, yextent);
3061 }
3062
3063 /* ----------- End of grid helpers ------------- */
3064
3065 /* vim: set shiftwidth=4 tabstop=8: */