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 (abs(eq[0]) < abs(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, 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, 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 * 2 * vec_x + vec_x;
1529     *yextent = height * vec_y;
1530 }
1531
1532 /* Doesn't use the previous method of generation, it pre-dates it!
1533  * A triangular grid is just about simple enough to do by "brute force" */
1534 static grid *grid_new_triangular(int width, int height, char *desc)
1535 {
1536     int x,y;
1537
1538     /* Vector for side of triangle - ratio is close to sqrt(3) */
1539     int vec_x = TRIANGLE_VEC_X;
1540     int vec_y = TRIANGLE_VEC_Y;
1541
1542     int index;
1543
1544     /* convenient alias */
1545     int w = width + 1;
1546
1547     grid *g = grid_empty();
1548     g->tilesize = TRIANGLE_TILESIZE;
1549
1550     g->num_faces = width * height * 2;
1551     g->num_dots = (width + 1) * (height + 1);
1552     g->faces = snewn(g->num_faces, grid_face);
1553     g->dots = snewn(g->num_dots, grid_dot);
1554
1555     /* generate dots */
1556     index = 0;
1557     for (y = 0; y <= height; y++) {
1558         for (x = 0; x <= width; x++) {
1559             grid_dot *d = g->dots + index;
1560             /* odd rows are offset to the right */
1561             d->order = 0;
1562             d->edges = NULL;
1563             d->faces = NULL;
1564             d->x = x * 2 * vec_x + ((y % 2) ? vec_x : 0);
1565             d->y = y * vec_y;
1566             index++;
1567         }
1568     }
1569
1570     /* generate faces */
1571     index = 0;
1572     for (y = 0; y < height; y++) {
1573         for (x = 0; x < width; x++) {
1574             /* initialise two faces for this (x,y) */
1575             grid_face *f1 = g->faces + index;
1576             grid_face *f2 = f1 + 1;
1577             f1->edges = NULL;
1578             f1->order = 3;
1579             f1->dots = snewn(f1->order, grid_dot*);
1580             f1->has_incentre = FALSE;
1581             f2->edges = NULL;
1582             f2->order = 3;
1583             f2->dots = snewn(f2->order, grid_dot*);
1584             f2->has_incentre = FALSE;
1585
1586             /* face descriptions depend on whether the row-number is
1587              * odd or even */
1588             if (y % 2) {
1589                 f1->dots[0] = g->dots + y       * w + x;
1590                 f1->dots[1] = g->dots + (y + 1) * w + x + 1;
1591                 f1->dots[2] = g->dots + (y + 1) * w + x;
1592                 f2->dots[0] = g->dots + y       * w + x;
1593                 f2->dots[1] = g->dots + y       * w + x + 1;
1594                 f2->dots[2] = g->dots + (y + 1) * w + x + 1;
1595             } else {
1596                 f1->dots[0] = g->dots + y       * w + x;
1597                 f1->dots[1] = g->dots + y       * w + x + 1;
1598                 f1->dots[2] = g->dots + (y + 1) * w + x;
1599                 f2->dots[0] = g->dots + y       * w + x + 1;
1600                 f2->dots[1] = g->dots + (y + 1) * w + x + 1;
1601                 f2->dots[2] = g->dots + (y + 1) * w + x;
1602             }
1603             index += 2;
1604         }
1605     }
1606
1607     grid_make_consistent(g);
1608     return g;
1609 }
1610
1611 #define SNUBSQUARE_TILESIZE 18
1612 /* Vector for side of triangle - ratio is close to sqrt(3) */
1613 #define SNUBSQUARE_A 15
1614 #define SNUBSQUARE_B 26
1615
1616 static void grid_size_snubsquare(int width, int height,
1617                           int *tilesize, int *xextent, int *yextent)
1618 {
1619     int a = SNUBSQUARE_A;
1620     int b = SNUBSQUARE_B;
1621
1622     *tilesize = SNUBSQUARE_TILESIZE;
1623     *xextent = (a+b) * (width-1) + a + b;
1624     *yextent = (a+b) * (height-1) + a + b;
1625 }
1626
1627 static grid *grid_new_snubsquare(int width, int height, char *desc)
1628 {
1629     int x, y;
1630     int a = SNUBSQUARE_A;
1631     int b = SNUBSQUARE_B;
1632
1633     /* Upper bounds - don't have to be exact */
1634     int max_faces = 3 * width * height;
1635     int max_dots = 2 * (width + 1) * (height + 1);
1636
1637     tree234 *points;
1638
1639     grid *g = grid_empty();
1640     g->tilesize = SNUBSQUARE_TILESIZE;
1641     g->faces = snewn(max_faces, grid_face);
1642     g->dots = snewn(max_dots, grid_dot);
1643
1644     points = newtree234(grid_point_cmp_fn);
1645
1646     for (y = 0; y < height; y++) {
1647         for (x = 0; x < width; x++) {
1648             grid_dot *d;
1649             /* face position */
1650             int px = (a + b) * x;
1651             int py = (a + b) * y;
1652
1653             /* generate square faces */
1654             grid_face_add_new(g, 4);
1655             if ((x + y) % 2) {
1656                 d = grid_get_dot(g, points, px + a, py);
1657                 grid_face_set_dot(g, d, 0);
1658                 d = grid_get_dot(g, points, px + a + b, py + a);
1659                 grid_face_set_dot(g, d, 1);
1660                 d = grid_get_dot(g, points, px + b, py + a + b);
1661                 grid_face_set_dot(g, d, 2);
1662                 d = grid_get_dot(g, points, px, py + b);
1663                 grid_face_set_dot(g, d, 3);
1664             } else {
1665                 d = grid_get_dot(g, points, px + b, py);
1666                 grid_face_set_dot(g, d, 0);
1667                 d = grid_get_dot(g, points, px + a + b, py + b);
1668                 grid_face_set_dot(g, d, 1);
1669                 d = grid_get_dot(g, points, px + a, py + a + b);
1670                 grid_face_set_dot(g, d, 2);
1671                 d = grid_get_dot(g, points, px, py + a);
1672                 grid_face_set_dot(g, d, 3);
1673             }
1674
1675             /* generate up/down triangles */
1676             if (x > 0) {
1677                 grid_face_add_new(g, 3);
1678                 if ((x + y) % 2) {
1679                     d = grid_get_dot(g, points, px + a, py);
1680                     grid_face_set_dot(g, d, 0);
1681                     d = grid_get_dot(g, points, px, py + b);
1682                     grid_face_set_dot(g, d, 1);
1683                     d = grid_get_dot(g, points, px - a, py);
1684                     grid_face_set_dot(g, d, 2);
1685                 } else {
1686                     d = grid_get_dot(g, points, px, py + a);
1687                     grid_face_set_dot(g, d, 0);
1688                     d = grid_get_dot(g, points, px + a, py + a + b);
1689                     grid_face_set_dot(g, d, 1);
1690                     d = grid_get_dot(g, points, px - a, py + a + b);
1691                     grid_face_set_dot(g, d, 2);
1692                 }
1693             }
1694
1695             /* generate left/right triangles */
1696             if (y > 0) {
1697                 grid_face_add_new(g, 3);
1698                 if ((x + y) % 2) {
1699                     d = grid_get_dot(g, points, px + a, py);
1700                     grid_face_set_dot(g, d, 0);
1701                     d = grid_get_dot(g, points, px + a + b, py - a);
1702                     grid_face_set_dot(g, d, 1);
1703                     d = grid_get_dot(g, points, px + a + b, py + a);
1704                     grid_face_set_dot(g, d, 2);
1705                 } else {
1706                     d = grid_get_dot(g, points, px, py - a);
1707                     grid_face_set_dot(g, d, 0);
1708                     d = grid_get_dot(g, points, px + b, py);
1709                     grid_face_set_dot(g, d, 1);
1710                     d = grid_get_dot(g, points, px, py + a);
1711                     grid_face_set_dot(g, d, 2);
1712                 }
1713             }
1714         }
1715     }
1716
1717     freetree234(points);
1718     assert(g->num_faces <= max_faces);
1719     assert(g->num_dots <= max_dots);
1720
1721     grid_make_consistent(g);
1722     return g;
1723 }
1724
1725 #define CAIRO_TILESIZE 40
1726 /* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
1727 #define CAIRO_A 14
1728 #define CAIRO_B 31
1729
1730 static void grid_size_cairo(int width, int height,
1731                           int *tilesize, int *xextent, int *yextent)
1732 {
1733     int b = CAIRO_B; /* a unused in determining grid size. */
1734
1735     *tilesize = CAIRO_TILESIZE;
1736     *xextent = 2*b*(width-1) + 2*b;
1737     *yextent = 2*b*(height-1) + 2*b;
1738 }
1739
1740 static grid *grid_new_cairo(int width, int height, char *desc)
1741 {
1742     int x, y;
1743     int a = CAIRO_A;
1744     int b = CAIRO_B;
1745
1746     /* Upper bounds - don't have to be exact */
1747     int max_faces = 2 * width * height;
1748     int max_dots = 3 * (width + 1) * (height + 1);
1749
1750     tree234 *points;
1751
1752     grid *g = grid_empty();
1753     g->tilesize = CAIRO_TILESIZE;
1754     g->faces = snewn(max_faces, grid_face);
1755     g->dots = snewn(max_dots, grid_dot);
1756
1757     points = newtree234(grid_point_cmp_fn);
1758
1759     for (y = 0; y < height; y++) {
1760         for (x = 0; x < width; x++) {
1761             grid_dot *d;
1762             /* cell position */
1763             int px = 2 * b * x;
1764             int py = 2 * b * y;
1765
1766             /* horizontal pentagons */
1767             if (y > 0) {
1768                 grid_face_add_new(g, 5);
1769                 if ((x + y) % 2) {
1770                     d = grid_get_dot(g, points, px + a, py - b);
1771                     grid_face_set_dot(g, d, 0);
1772                     d = grid_get_dot(g, points, px + 2*b - a, py - b);
1773                     grid_face_set_dot(g, d, 1);
1774                     d = grid_get_dot(g, points, px + 2*b, py);
1775                     grid_face_set_dot(g, d, 2);
1776                     d = grid_get_dot(g, points, px + b, py + a);
1777                     grid_face_set_dot(g, d, 3);
1778                     d = grid_get_dot(g, points, px, py);
1779                     grid_face_set_dot(g, d, 4);
1780                 } else {
1781                     d = grid_get_dot(g, points, px, py);
1782                     grid_face_set_dot(g, d, 0);
1783                     d = grid_get_dot(g, points, px + b, py - a);
1784                     grid_face_set_dot(g, d, 1);
1785                     d = grid_get_dot(g, points, px + 2*b, py);
1786                     grid_face_set_dot(g, d, 2);
1787                     d = grid_get_dot(g, points, px + 2*b - a, py + b);
1788                     grid_face_set_dot(g, d, 3);
1789                     d = grid_get_dot(g, points, px + a, py + b);
1790                     grid_face_set_dot(g, d, 4);
1791                 }
1792             }
1793             /* vertical pentagons */
1794             if (x > 0) {
1795                 grid_face_add_new(g, 5);
1796                 if ((x + y) % 2) {
1797                     d = grid_get_dot(g, points, px, py);
1798                     grid_face_set_dot(g, d, 0);
1799                     d = grid_get_dot(g, points, px + b, py + a);
1800                     grid_face_set_dot(g, d, 1);
1801                     d = grid_get_dot(g, points, px + b, py + 2*b - a);
1802                     grid_face_set_dot(g, d, 2);
1803                     d = grid_get_dot(g, points, px, py + 2*b);
1804                     grid_face_set_dot(g, d, 3);
1805                     d = grid_get_dot(g, points, px - a, py + b);
1806                     grid_face_set_dot(g, d, 4);
1807                 } else {
1808                     d = grid_get_dot(g, points, px, py);
1809                     grid_face_set_dot(g, d, 0);
1810                     d = grid_get_dot(g, points, px + a, py + b);
1811                     grid_face_set_dot(g, d, 1);
1812                     d = grid_get_dot(g, points, px, py + 2*b);
1813                     grid_face_set_dot(g, d, 2);
1814                     d = grid_get_dot(g, points, px - b, py + 2*b - a);
1815                     grid_face_set_dot(g, d, 3);
1816                     d = grid_get_dot(g, points, px - b, py + a);
1817                     grid_face_set_dot(g, d, 4);
1818                 }
1819             }
1820         }
1821     }
1822
1823     freetree234(points);
1824     assert(g->num_faces <= max_faces);
1825     assert(g->num_dots <= max_dots);
1826
1827     grid_make_consistent(g);
1828     return g;
1829 }
1830
1831 #define GREATHEX_TILESIZE 18
1832 /* Vector for side of triangle - ratio is close to sqrt(3) */
1833 #define GREATHEX_A 15
1834 #define GREATHEX_B 26
1835
1836 static void grid_size_greathexagonal(int width, int height,
1837                           int *tilesize, int *xextent, int *yextent)
1838 {
1839     int a = GREATHEX_A;
1840     int b = GREATHEX_B;
1841
1842     *tilesize = GREATHEX_TILESIZE;
1843     *xextent = (3*a + b) * (width-1) + 4*a;
1844     *yextent = (2*a + 2*b) * (height-1) + 3*b + a;
1845 }
1846
1847 static grid *grid_new_greathexagonal(int width, int height, char *desc)
1848 {
1849     int x, y;
1850     int a = GREATHEX_A;
1851     int b = GREATHEX_B;
1852
1853     /* Upper bounds - don't have to be exact */
1854     int max_faces = 6 * (width + 1) * (height + 1);
1855     int max_dots = 6 * width * height;
1856
1857     tree234 *points;
1858
1859     grid *g = grid_empty();
1860     g->tilesize = GREATHEX_TILESIZE;
1861     g->faces = snewn(max_faces, grid_face);
1862     g->dots = snewn(max_dots, grid_dot);
1863
1864     points = newtree234(grid_point_cmp_fn);
1865
1866     for (y = 0; y < height; y++) {
1867         for (x = 0; x < width; x++) {
1868             grid_dot *d;
1869             /* centre of hexagon */
1870             int px = (3*a + b) * x;
1871             int py = (2*a + 2*b) * y;
1872             if (x % 2)
1873                 py += a + b;
1874
1875             /* hexagon */
1876             grid_face_add_new(g, 6);
1877             d = grid_get_dot(g, points, px - a, py - b);
1878             grid_face_set_dot(g, d, 0);
1879             d = grid_get_dot(g, points, px + a, py - b);
1880             grid_face_set_dot(g, d, 1);
1881             d = grid_get_dot(g, points, px + 2*a, py);
1882             grid_face_set_dot(g, d, 2);
1883             d = grid_get_dot(g, points, px + a, py + b);
1884             grid_face_set_dot(g, d, 3);
1885             d = grid_get_dot(g, points, px - a, py + b);
1886             grid_face_set_dot(g, d, 4);
1887             d = grid_get_dot(g, points, px - 2*a, py);
1888             grid_face_set_dot(g, d, 5);
1889
1890             /* square below hexagon */
1891             if (y < height - 1) {
1892                 grid_face_add_new(g, 4);
1893                 d = grid_get_dot(g, points, px - a, py + b);
1894                 grid_face_set_dot(g, d, 0);
1895                 d = grid_get_dot(g, points, px + a, py + b);
1896                 grid_face_set_dot(g, d, 1);
1897                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
1898                 grid_face_set_dot(g, d, 2);
1899                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
1900                 grid_face_set_dot(g, d, 3);
1901             }
1902
1903             /* square below right */
1904             if ((x < width - 1) && (((x % 2) == 0) || (y < height - 1))) {
1905                 grid_face_add_new(g, 4);
1906                 d = grid_get_dot(g, points, px + 2*a, py);
1907                 grid_face_set_dot(g, d, 0);
1908                 d = grid_get_dot(g, points, px + 2*a + b, py + a);
1909                 grid_face_set_dot(g, d, 1);
1910                 d = grid_get_dot(g, points, px + a + b, py + a + b);
1911                 grid_face_set_dot(g, d, 2);
1912                 d = grid_get_dot(g, points, px + a, py + b);
1913                 grid_face_set_dot(g, d, 3);
1914             }
1915
1916             /* square below left */
1917             if ((x > 0) && (((x % 2) == 0) || (y < height - 1))) {
1918                 grid_face_add_new(g, 4);
1919                 d = grid_get_dot(g, points, px - 2*a, py);
1920                 grid_face_set_dot(g, d, 0);
1921                 d = grid_get_dot(g, points, px - a, py + b);
1922                 grid_face_set_dot(g, d, 1);
1923                 d = grid_get_dot(g, points, px - a - b, py + a + b);
1924                 grid_face_set_dot(g, d, 2);
1925                 d = grid_get_dot(g, points, px - 2*a - b, py + a);
1926                 grid_face_set_dot(g, d, 3);
1927             }
1928
1929             /* Triangle below right */
1930             if ((x < width - 1) && (y < height - 1)) {
1931                 grid_face_add_new(g, 3);
1932                 d = grid_get_dot(g, points, px + a, py + b);
1933                 grid_face_set_dot(g, d, 0);
1934                 d = grid_get_dot(g, points, px + a + b, py + a + b);
1935                 grid_face_set_dot(g, d, 1);
1936                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
1937                 grid_face_set_dot(g, d, 2);
1938             }
1939
1940             /* Triangle below left */
1941             if ((x > 0) && (y < height - 1)) {
1942                 grid_face_add_new(g, 3);
1943                 d = grid_get_dot(g, points, px - a, py + b);
1944                 grid_face_set_dot(g, d, 0);
1945                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
1946                 grid_face_set_dot(g, d, 1);
1947                 d = grid_get_dot(g, points, px - a - b, py + a + b);
1948                 grid_face_set_dot(g, d, 2);
1949             }
1950         }
1951     }
1952
1953     freetree234(points);
1954     assert(g->num_faces <= max_faces);
1955     assert(g->num_dots <= max_dots);
1956
1957     grid_make_consistent(g);
1958     return g;
1959 }
1960
1961 #define OCTAGONAL_TILESIZE 40
1962 /* b/a approx sqrt(2) */
1963 #define OCTAGONAL_A 29
1964 #define OCTAGONAL_B 41
1965
1966 static void grid_size_octagonal(int width, int height,
1967                           int *tilesize, int *xextent, int *yextent)
1968 {
1969     int a = OCTAGONAL_A;
1970     int b = OCTAGONAL_B;
1971
1972     *tilesize = OCTAGONAL_TILESIZE;
1973     *xextent = (2*a + b) * width;
1974     *yextent = (2*a + b) * height;
1975 }
1976
1977 static grid *grid_new_octagonal(int width, int height, char *desc)
1978 {
1979     int x, y;
1980     int a = OCTAGONAL_A;
1981     int b = OCTAGONAL_B;
1982
1983     /* Upper bounds - don't have to be exact */
1984     int max_faces = 2 * width * height;
1985     int max_dots = 4 * (width + 1) * (height + 1);
1986
1987     tree234 *points;
1988
1989     grid *g = grid_empty();
1990     g->tilesize = OCTAGONAL_TILESIZE;
1991     g->faces = snewn(max_faces, grid_face);
1992     g->dots = snewn(max_dots, grid_dot);
1993
1994     points = newtree234(grid_point_cmp_fn);
1995
1996     for (y = 0; y < height; y++) {
1997         for (x = 0; x < width; x++) {
1998             grid_dot *d;
1999             /* cell position */
2000             int px = (2*a + b) * x;
2001             int py = (2*a + b) * y;
2002             /* octagon */
2003             grid_face_add_new(g, 8);
2004             d = grid_get_dot(g, points, px + a, py);
2005             grid_face_set_dot(g, d, 0);
2006             d = grid_get_dot(g, points, px + a + b, py);
2007             grid_face_set_dot(g, d, 1);
2008             d = grid_get_dot(g, points, px + 2*a + b, py + a);
2009             grid_face_set_dot(g, d, 2);
2010             d = grid_get_dot(g, points, px + 2*a + b, py + a + b);
2011             grid_face_set_dot(g, d, 3);
2012             d = grid_get_dot(g, points, px + a + b, py + 2*a + b);
2013             grid_face_set_dot(g, d, 4);
2014             d = grid_get_dot(g, points, px + a, py + 2*a + b);
2015             grid_face_set_dot(g, d, 5);
2016             d = grid_get_dot(g, points, px, py + a + b);
2017             grid_face_set_dot(g, d, 6);
2018             d = grid_get_dot(g, points, px, py + a);
2019             grid_face_set_dot(g, d, 7);
2020
2021             /* diamond */
2022             if ((x > 0) && (y > 0)) {
2023                 grid_face_add_new(g, 4);
2024                 d = grid_get_dot(g, points, px, py - a);
2025                 grid_face_set_dot(g, d, 0);
2026                 d = grid_get_dot(g, points, px + a, py);
2027                 grid_face_set_dot(g, d, 1);
2028                 d = grid_get_dot(g, points, px, py + a);
2029                 grid_face_set_dot(g, d, 2);
2030                 d = grid_get_dot(g, points, px - a, py);
2031                 grid_face_set_dot(g, d, 3);
2032             }
2033         }
2034     }
2035
2036     freetree234(points);
2037     assert(g->num_faces <= max_faces);
2038     assert(g->num_dots <= max_dots);
2039
2040     grid_make_consistent(g);
2041     return g;
2042 }
2043
2044 #define KITE_TILESIZE 40
2045 /* b/a approx sqrt(3) */
2046 #define KITE_A 15
2047 #define KITE_B 26
2048
2049 static void grid_size_kites(int width, int height,
2050                      int *tilesize, int *xextent, int *yextent)
2051 {
2052     int a = KITE_A;
2053     int b = KITE_B;
2054
2055     *tilesize = KITE_TILESIZE;
2056     *xextent = 4*b * width + 2*b;
2057     *yextent = 6*a * (height-1) + 8*a;
2058 }
2059
2060 static grid *grid_new_kites(int width, int height, char *desc)
2061 {
2062     int x, y;
2063     int a = KITE_A;
2064     int b = KITE_B;
2065
2066     /* Upper bounds - don't have to be exact */
2067     int max_faces = 6 * width * height;
2068     int max_dots = 6 * (width + 1) * (height + 1);
2069
2070     tree234 *points;
2071
2072     grid *g = grid_empty();
2073     g->tilesize = KITE_TILESIZE;
2074     g->faces = snewn(max_faces, grid_face);
2075     g->dots = snewn(max_dots, grid_dot);
2076
2077     points = newtree234(grid_point_cmp_fn);
2078
2079     for (y = 0; y < height; y++) {
2080         for (x = 0; x < width; x++) {
2081             grid_dot *d;
2082             /* position of order-6 dot */
2083             int px = 4*b * x;
2084             int py = 6*a * y;
2085             if (y % 2)
2086                 px += 2*b;
2087
2088             /* kite pointing up-left */
2089             grid_face_add_new(g, 4);
2090             d = grid_get_dot(g, points, px, py);
2091             grid_face_set_dot(g, d, 0);
2092             d = grid_get_dot(g, points, px + 2*b, py);
2093             grid_face_set_dot(g, d, 1);
2094             d = grid_get_dot(g, points, px + 2*b, py + 2*a);
2095             grid_face_set_dot(g, d, 2);
2096             d = grid_get_dot(g, points, px + b, py + 3*a);
2097             grid_face_set_dot(g, d, 3);
2098
2099             /* kite pointing up */
2100             grid_face_add_new(g, 4);
2101             d = grid_get_dot(g, points, px, py);
2102             grid_face_set_dot(g, d, 0);
2103             d = grid_get_dot(g, points, px + b, py + 3*a);
2104             grid_face_set_dot(g, d, 1);
2105             d = grid_get_dot(g, points, px, py + 4*a);
2106             grid_face_set_dot(g, d, 2);
2107             d = grid_get_dot(g, points, px - b, py + 3*a);
2108             grid_face_set_dot(g, d, 3);
2109
2110             /* kite pointing up-right */
2111             grid_face_add_new(g, 4);
2112             d = grid_get_dot(g, points, px, py);
2113             grid_face_set_dot(g, d, 0);
2114             d = grid_get_dot(g, points, px - b, py + 3*a);
2115             grid_face_set_dot(g, d, 1);
2116             d = grid_get_dot(g, points, px - 2*b, py + 2*a);
2117             grid_face_set_dot(g, d, 2);
2118             d = grid_get_dot(g, points, px - 2*b, py);
2119             grid_face_set_dot(g, d, 3);
2120
2121             /* kite pointing down-right */
2122             grid_face_add_new(g, 4);
2123             d = grid_get_dot(g, points, px, py);
2124             grid_face_set_dot(g, d, 0);
2125             d = grid_get_dot(g, points, px - 2*b, py);
2126             grid_face_set_dot(g, d, 1);
2127             d = grid_get_dot(g, points, px - 2*b, py - 2*a);
2128             grid_face_set_dot(g, d, 2);
2129             d = grid_get_dot(g, points, px - b, py - 3*a);
2130             grid_face_set_dot(g, d, 3);
2131
2132             /* kite pointing down */
2133             grid_face_add_new(g, 4);
2134             d = grid_get_dot(g, points, px, py);
2135             grid_face_set_dot(g, d, 0);
2136             d = grid_get_dot(g, points, px - b, py - 3*a);
2137             grid_face_set_dot(g, d, 1);
2138             d = grid_get_dot(g, points, px, py - 4*a);
2139             grid_face_set_dot(g, d, 2);
2140             d = grid_get_dot(g, points, px + b, py - 3*a);
2141             grid_face_set_dot(g, d, 3);
2142
2143             /* kite pointing down-left */
2144             grid_face_add_new(g, 4);
2145             d = grid_get_dot(g, points, px, py);
2146             grid_face_set_dot(g, d, 0);
2147             d = grid_get_dot(g, points, px + b, py - 3*a);
2148             grid_face_set_dot(g, d, 1);
2149             d = grid_get_dot(g, points, px + 2*b, py - 2*a);
2150             grid_face_set_dot(g, d, 2);
2151             d = grid_get_dot(g, points, px + 2*b, py);
2152             grid_face_set_dot(g, d, 3);
2153         }
2154     }
2155
2156     freetree234(points);
2157     assert(g->num_faces <= max_faces);
2158     assert(g->num_dots <= max_dots);
2159
2160     grid_make_consistent(g);
2161     return g;
2162 }
2163
2164 #define FLORET_TILESIZE 150
2165 /* -py/px is close to tan(30 - atan(sqrt(3)/9))
2166  * using py=26 makes everything lean to the left, rather than right
2167  */
2168 #define FLORET_PX 75
2169 #define FLORET_PY -26
2170
2171 static void grid_size_floret(int width, int height,
2172                           int *tilesize, int *xextent, int *yextent)
2173 {
2174     int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
2175     int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
2176     int ry = qy-py;
2177     /* rx unused in determining grid size. */
2178
2179     *tilesize = FLORET_TILESIZE;
2180     *xextent = (6*px+3*qx)/2 * (width-1) + 4*qx + 2*px;
2181     *yextent = (5*qy-4*py) * (height-1) + 4*qy + 2*ry;
2182 }
2183
2184 static grid *grid_new_floret(int width, int height, char *desc)
2185 {
2186     int x, y;
2187     /* Vectors for sides; weird numbers needed to keep puzzle aligned with window
2188      * -py/px is close to tan(30 - atan(sqrt(3)/9))
2189      * using py=26 makes everything lean to the left, rather than right
2190      */
2191     int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
2192     int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
2193     int rx = qx-px, ry = qy-py;                 /* |(-15,  78)| = 79.38 */
2194
2195     /* Upper bounds - don't have to be exact */
2196     int max_faces = 6 * width * height;
2197     int max_dots = 9 * (width + 1) * (height + 1);
2198
2199     tree234 *points;
2200
2201     grid *g = grid_empty();
2202     g->tilesize = FLORET_TILESIZE;
2203     g->faces = snewn(max_faces, grid_face);
2204     g->dots = snewn(max_dots, grid_dot);
2205
2206     points = newtree234(grid_point_cmp_fn);
2207
2208     /* generate pentagonal faces */
2209     for (y = 0; y < height; y++) {
2210         for (x = 0; x < width; x++) {
2211             grid_dot *d;
2212             /* face centre */
2213             int cx = (6*px+3*qx)/2 * x;
2214             int cy = (4*py-5*qy) * y;
2215             if (x % 2)
2216                 cy -= (4*py-5*qy)/2;
2217             else if (y && y == height-1)
2218                 continue; /* make better looking grids?  try 3x3 for instance */
2219
2220             grid_face_add_new(g, 5);
2221             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2222             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 1);
2223             d = grid_get_dot(g, points, cx+2*rx+qx, cy+2*ry+qy); grid_face_set_dot(g, d, 2);
2224             d = grid_get_dot(g, points, cx+2*qx+rx, cy+2*qy+ry); grid_face_set_dot(g, d, 3);
2225             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 4);
2226
2227             grid_face_add_new(g, 5);
2228             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2229             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 1);
2230             d = grid_get_dot(g, points, cx+2*qx+px, cy+2*qy+py); grid_face_set_dot(g, d, 2);
2231             d = grid_get_dot(g, points, cx+2*px+qx, cy+2*py+qy); grid_face_set_dot(g, d, 3);
2232             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 4);
2233
2234             grid_face_add_new(g, 5);
2235             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2236             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 1);
2237             d = grid_get_dot(g, points, cx+2*px-rx, cy+2*py-ry); grid_face_set_dot(g, d, 2);
2238             d = grid_get_dot(g, points, cx-2*rx+px, cy-2*ry+py); grid_face_set_dot(g, d, 3);
2239             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 4);
2240
2241             grid_face_add_new(g, 5);
2242             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2243             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 1);
2244             d = grid_get_dot(g, points, cx-2*rx-qx, cy-2*ry-qy); grid_face_set_dot(g, d, 2);
2245             d = grid_get_dot(g, points, cx-2*qx-rx, cy-2*qy-ry); grid_face_set_dot(g, d, 3);
2246             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 4);
2247
2248             grid_face_add_new(g, 5);
2249             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2250             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 1);
2251             d = grid_get_dot(g, points, cx-2*qx-px, cy-2*qy-py); grid_face_set_dot(g, d, 2);
2252             d = grid_get_dot(g, points, cx-2*px-qx, cy-2*py-qy); grid_face_set_dot(g, d, 3);
2253             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 4);
2254
2255             grid_face_add_new(g, 5);
2256             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
2257             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 1);
2258             d = grid_get_dot(g, points, cx-2*px+rx, cy-2*py+ry); grid_face_set_dot(g, d, 2);
2259             d = grid_get_dot(g, points, cx+2*rx-px, cy+2*ry-py); grid_face_set_dot(g, d, 3);
2260             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 4);
2261         }
2262     }
2263
2264     freetree234(points);
2265     assert(g->num_faces <= max_faces);
2266     assert(g->num_dots <= max_dots);
2267
2268     grid_make_consistent(g);
2269     return g;
2270 }
2271
2272 /* DODEC_* are used for dodecagonal and great-dodecagonal grids. */
2273 #define DODEC_TILESIZE 26
2274 /* Vector for side of triangle - ratio is close to sqrt(3) */
2275 #define DODEC_A 15
2276 #define DODEC_B 26
2277
2278 static void grid_size_dodecagonal(int width, int height,
2279                           int *tilesize, int *xextent, int *yextent)
2280 {
2281     int a = DODEC_A;
2282     int b = DODEC_B;
2283
2284     *tilesize = DODEC_TILESIZE;
2285     *xextent = (4*a + 2*b) * (width-1) + 3*(2*a + b);
2286     *yextent = (3*a + 2*b) * (height-1) + 2*(2*a + b);
2287 }
2288
2289 static grid *grid_new_dodecagonal(int width, int height, char *desc)
2290 {
2291     int x, y;
2292     int a = DODEC_A;
2293     int b = DODEC_B;
2294
2295     /* Upper bounds - don't have to be exact */
2296     int max_faces = 3 * width * height;
2297     int max_dots = 14 * width * height;
2298
2299     tree234 *points;
2300
2301     grid *g = grid_empty();
2302     g->tilesize = DODEC_TILESIZE;
2303     g->faces = snewn(max_faces, grid_face);
2304     g->dots = snewn(max_dots, grid_dot);
2305
2306     points = newtree234(grid_point_cmp_fn);
2307
2308     for (y = 0; y < height; y++) {
2309         for (x = 0; x < width; x++) {
2310             grid_dot *d;
2311             /* centre of dodecagon */
2312             int px = (4*a + 2*b) * x;
2313             int py = (3*a + 2*b) * y;
2314             if (y % 2)
2315                 px += 2*a + b;
2316
2317             /* dodecagon */
2318             grid_face_add_new(g, 12);
2319             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
2320             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
2321             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
2322             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
2323             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
2324             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
2325             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
2326             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
2327             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
2328             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
2329             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
2330             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
2331
2332             /* triangle below dodecagon */
2333             if ((y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
2334                 grid_face_add_new(g, 3);
2335                 d = grid_get_dot(g, points, px + a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2336                 d = grid_get_dot(g, points, px    , py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2337                 d = grid_get_dot(g, points, px - a, py + (2*a +   b)); grid_face_set_dot(g, d, 2);
2338             }
2339
2340             /* triangle above dodecagon */
2341             if ((y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
2342                 grid_face_add_new(g, 3);
2343                 d = grid_get_dot(g, points, px - a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2344                 d = grid_get_dot(g, points, px    , py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2345                 d = grid_get_dot(g, points, px + a, py - (2*a +   b)); grid_face_set_dot(g, d, 2);
2346             }
2347         }
2348     }
2349
2350     freetree234(points);
2351     assert(g->num_faces <= max_faces);
2352     assert(g->num_dots <= max_dots);
2353
2354     grid_make_consistent(g);
2355     return g;
2356 }
2357
2358 static void grid_size_greatdodecagonal(int width, int height,
2359                           int *tilesize, int *xextent, int *yextent)
2360 {
2361     int a = DODEC_A;
2362     int b = DODEC_B;
2363
2364     *tilesize = DODEC_TILESIZE;
2365     *xextent = (6*a + 2*b) * (width-1) + 2*(2*a + b) + 3*a + b;
2366     *yextent = (3*a + 3*b) * (height-1) + 2*(2*a + b);
2367 }
2368
2369 static grid *grid_new_greatdodecagonal(int width, int height, char *desc)
2370 {
2371     int x, y;
2372     /* Vector for side of triangle - ratio is close to sqrt(3) */
2373     int a = DODEC_A;
2374     int b = DODEC_B;
2375
2376     /* Upper bounds - don't have to be exact */
2377     int max_faces = 30 * width * height;
2378     int max_dots = 200 * width * height;
2379
2380     tree234 *points;
2381
2382     grid *g = grid_empty();
2383     g->tilesize = DODEC_TILESIZE;
2384     g->faces = snewn(max_faces, grid_face);
2385     g->dots = snewn(max_dots, grid_dot);
2386
2387     points = newtree234(grid_point_cmp_fn);
2388
2389     for (y = 0; y < height; y++) {
2390         for (x = 0; x < width; x++) {
2391             grid_dot *d;
2392             /* centre of dodecagon */
2393             int px = (6*a + 2*b) * x;
2394             int py = (3*a + 3*b) * y;
2395             if (y % 2)
2396                 px += 3*a + b;
2397
2398             /* dodecagon */
2399             grid_face_add_new(g, 12);
2400             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
2401             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
2402             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
2403             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
2404             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
2405             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
2406             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
2407             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
2408             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
2409             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
2410             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
2411             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
2412
2413             /* hexagon below dodecagon */
2414             if (y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2415                 grid_face_add_new(g, 6);
2416                 d = grid_get_dot(g, points, px +   a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2417                 d = grid_get_dot(g, points, px + 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2418                 d = grid_get_dot(g, points, px +   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2419                 d = grid_get_dot(g, points, px -   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2420                 d = grid_get_dot(g, points, px - 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2421                 d = grid_get_dot(g, points, px -   a, py + (2*a +   b)); grid_face_set_dot(g, d, 5);
2422             }
2423
2424             /* hexagon above dodecagon */
2425             if (y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2426                 grid_face_add_new(g, 6);
2427                 d = grid_get_dot(g, points, px -   a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2428                 d = grid_get_dot(g, points, px - 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2429                 d = grid_get_dot(g, points, px -   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2430                 d = grid_get_dot(g, points, px +   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2431                 d = grid_get_dot(g, points, px + 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2432                 d = grid_get_dot(g, points, px +   a, py - (2*a +   b)); grid_face_set_dot(g, d, 5);
2433             }
2434
2435             /* square on right of dodecagon */
2436             if (x < width - 1) {
2437                 grid_face_add_new(g, 4);
2438                 d = grid_get_dot(g, points, px + 2*a + b, py - a); grid_face_set_dot(g, d, 0);
2439                 d = grid_get_dot(g, points, px + 4*a + b, py - a); grid_face_set_dot(g, d, 1);
2440                 d = grid_get_dot(g, points, px + 4*a + b, py + a); grid_face_set_dot(g, d, 2);
2441                 d = grid_get_dot(g, points, px + 2*a + b, py + a); grid_face_set_dot(g, d, 3);
2442             }
2443
2444             /* square on top right of dodecagon */
2445             if (y && (x < width - 1 || !(y % 2))) {
2446                 grid_face_add_new(g, 4);
2447                 d = grid_get_dot(g, points, px + (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2448                 d = grid_get_dot(g, points, px + (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2449                 d = grid_get_dot(g, points, px + (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 2);
2450                 d = grid_get_dot(g, points, px + (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 3);
2451             }
2452
2453             /* square on top left of dodecagon */
2454             if (y && (x || (y % 2))) {
2455                 grid_face_add_new(g, 4);
2456                 d = grid_get_dot(g, points, px - (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 0);
2457                 d = grid_get_dot(g, points, px - (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 1);
2458                 d = grid_get_dot(g, points, px - (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 2);
2459                 d = grid_get_dot(g, points, px - (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 3);
2460             }
2461         }
2462     }
2463
2464     freetree234(points);
2465     assert(g->num_faces <= max_faces);
2466     assert(g->num_dots <= max_dots);
2467
2468     grid_make_consistent(g);
2469     return g;
2470 }
2471
2472 typedef struct setface_ctx
2473 {
2474     int xmin, xmax, ymin, ymax;
2475
2476     grid *g;
2477     tree234 *points;
2478 } setface_ctx;
2479
2480 static double round_int_nearest_away(double r)
2481 {
2482     return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
2483 }
2484
2485 static int set_faces(penrose_state *state, vector *vs, int n, int depth)
2486 {
2487     setface_ctx *sf_ctx = (setface_ctx *)state->ctx;
2488     int i;
2489     int xs[4], ys[4];
2490
2491     if (depth < state->max_depth) return 0;
2492 #ifdef DEBUG_PENROSE
2493     if (n != 4) return 0; /* triangles are sent as debugging. */
2494 #endif
2495
2496     for (i = 0; i < n; i++) {
2497         double tx = v_x(vs, i), ty = v_y(vs, i);
2498
2499         xs[i] = (int)round_int_nearest_away(tx);
2500         ys[i] = (int)round_int_nearest_away(ty);
2501
2502         if (xs[i] < sf_ctx->xmin || xs[i] > sf_ctx->xmax) return 0;
2503         if (ys[i] < sf_ctx->ymin || ys[i] > sf_ctx->ymax) return 0;
2504     }
2505
2506     grid_face_add_new(sf_ctx->g, n);
2507     debug(("penrose: new face l=%f gen=%d...",
2508            penrose_side_length(state->start_size, depth), depth));
2509     for (i = 0; i < n; i++) {
2510         grid_dot *d = grid_get_dot(sf_ctx->g, sf_ctx->points,
2511                                    xs[i], ys[i]);
2512         grid_face_set_dot(sf_ctx->g, d, i);
2513         debug((" ... dot 0x%x (%d,%d) (was %2.2f,%2.2f)",
2514                d, d->x, d->y, v_x(vs, i), v_y(vs, i)));
2515     }
2516
2517     return 0;
2518 }
2519
2520 #define PENROSE_TILESIZE 100
2521
2522 static void grid_size_penrose(int width, int height,
2523                        int *tilesize, int *xextent, int *yextent)
2524 {
2525     int l = PENROSE_TILESIZE;
2526
2527     *tilesize = l;
2528     *xextent = l * width;
2529     *yextent = l * height;
2530 }
2531
2532 static char *grid_new_desc_penrose(grid_type type, int width, int height, random_state *rs)
2533 {
2534     int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff;
2535     double outer_radius;
2536     int inner_radius;
2537     char gd[255];
2538     int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
2539
2540     /* We want to produce a random bit of penrose tiling, so we calculate
2541      * a random offset (within the patch that penrose.c calculates for us)
2542      * and an angle (multiple of 36) to rotate the patch. */
2543
2544     penrose_calculate_size(which, tilesize, width, height,
2545                            &outer_radius, &startsz, &depth);
2546
2547     /* Calculate radius of (circumcircle of) patch, subtract from
2548      * radius calculated. */
2549     inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
2550
2551     /* Pick a random offset (the easy way: choose within outer square,
2552      * discarding while it's outside the circle) */
2553     do {
2554         xoff = random_upto(rs, 2*inner_radius) - inner_radius;
2555         yoff = random_upto(rs, 2*inner_radius) - inner_radius;
2556     } while (sqrt(xoff*xoff+yoff*yoff) > inner_radius);
2557
2558     aoff = random_upto(rs, 360/36) * 36;
2559
2560     debug(("grid_desc: ts %d, %dx%d patch, orad %2.2f irad %d",
2561            tilesize, width, height, outer_radius, inner_radius));
2562     debug(("    -> xoff %d yoff %d aoff %d", xoff, yoff, aoff));
2563
2564     sprintf(gd, "G%d,%d,%d", xoff, yoff, aoff);
2565
2566     return dupstr(gd);
2567 }
2568
2569 static char *grid_validate_desc_penrose(grid_type type, int width, int height, char *desc)
2570 {
2571     int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff, inner_radius;
2572     double outer_radius;
2573     int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
2574
2575     if (!desc)
2576         return "Missing grid description string.";
2577
2578     penrose_calculate_size(which, tilesize, width, height,
2579                            &outer_radius, &startsz, &depth);
2580     inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
2581
2582     if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
2583         return "Invalid format grid description string.";
2584
2585     if (sqrt(xoff*xoff + yoff*yoff) > inner_radius)
2586         return "Patch offset out of bounds.";
2587     if ((aoff % 36) != 0 || aoff < 0 || aoff >= 360)
2588         return "Angle offset out of bounds.";
2589
2590     return NULL;
2591 }
2592
2593 /*
2594  * We're asked for a grid of a particular size, and we generate enough
2595  * of the tiling so we can be sure to have enough random grid from which
2596  * to pick.
2597  */
2598
2599 static grid *grid_new_penrose(int width, int height, int which, char *desc)
2600 {
2601     int max_faces, max_dots, tilesize = PENROSE_TILESIZE;
2602     int xsz, ysz, xoff, yoff, aoff;
2603     double rradius;
2604
2605     tree234 *points;
2606     grid *g;
2607
2608     penrose_state ps;
2609     setface_ctx sf_ctx;
2610
2611     penrose_calculate_size(which, tilesize, width, height,
2612                            &rradius, &ps.start_size, &ps.max_depth);
2613
2614     debug(("penrose: w%d h%d, tile size %d, start size %d, depth %d",
2615            width, height, tilesize, ps.start_size, ps.max_depth));
2616
2617     ps.new_tile = set_faces;
2618     ps.ctx = &sf_ctx;
2619
2620     max_faces = (width*3) * (height*3); /* somewhat paranoid... */
2621     max_dots = max_faces * 4; /* ditto... */
2622
2623     g = grid_empty();
2624     g->tilesize = tilesize;
2625     g->faces = snewn(max_faces, grid_face);
2626     g->dots = snewn(max_dots, grid_dot);
2627
2628     points = newtree234(grid_point_cmp_fn);
2629
2630     memset(&sf_ctx, 0, sizeof(sf_ctx));
2631     sf_ctx.g = g;
2632     sf_ctx.points = points;
2633
2634     if (desc != NULL) {
2635         if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
2636             assert(!"Invalid grid description.");
2637     } else {
2638         xoff = yoff = aoff = 0;
2639     }
2640
2641     xsz = width * tilesize;
2642     ysz = height * tilesize;
2643
2644     sf_ctx.xmin = xoff - xsz/2;
2645     sf_ctx.xmax = xoff + xsz/2;
2646     sf_ctx.ymin = yoff - ysz/2;
2647     sf_ctx.ymax = yoff + ysz/2;
2648
2649     debug(("penrose: centre (%f, %f) xsz %f ysz %f",
2650            0.0, 0.0, xsz, ysz));
2651     debug(("penrose: x range (%f --> %f), y range (%f --> %f)",
2652            sf_ctx.xmin, sf_ctx.xmax, sf_ctx.ymin, sf_ctx.ymax));
2653
2654     penrose(&ps, which, aoff);
2655
2656     freetree234(points);
2657     assert(g->num_faces <= max_faces);
2658     assert(g->num_dots <= max_dots);
2659
2660     debug(("penrose: %d faces total (equivalent to %d wide by %d high)",
2661            g->num_faces, g->num_faces/height, g->num_faces/width));
2662
2663     grid_trim_vigorously(g);
2664     grid_make_consistent(g);
2665
2666     /*
2667      * Centre the grid in its originally promised rectangle.
2668      */
2669     g->lowest_x -= ((sf_ctx.xmax - sf_ctx.xmin) -
2670                     (g->highest_x - g->lowest_x)) / 2;
2671     g->highest_x = g->lowest_x + (sf_ctx.xmax - sf_ctx.xmin);
2672     g->lowest_y -= ((sf_ctx.ymax - sf_ctx.ymin) -
2673                     (g->highest_y - g->lowest_y)) / 2;
2674     g->highest_y = g->lowest_y + (sf_ctx.ymax - sf_ctx.ymin);
2675
2676     return g;
2677 }
2678
2679 static void grid_size_penrose_p2_kite(int width, int height,
2680                        int *tilesize, int *xextent, int *yextent)
2681 {
2682     grid_size_penrose(width, height, tilesize, xextent, yextent);
2683 }
2684
2685 static void grid_size_penrose_p3_thick(int width, int height,
2686                        int *tilesize, int *xextent, int *yextent)
2687 {
2688     grid_size_penrose(width, height, tilesize, xextent, yextent);
2689 }
2690
2691 static grid *grid_new_penrose_p2_kite(int width, int height, char *desc)
2692 {
2693     return grid_new_penrose(width, height, PENROSE_P2, desc);
2694 }
2695
2696 static grid *grid_new_penrose_p3_thick(int width, int height, char *desc)
2697 {
2698     return grid_new_penrose(width, height, PENROSE_P3, desc);
2699 }
2700
2701 /* ----------- End of grid generators ------------- */
2702
2703 #define FNNEW(upper,lower) &grid_new_ ## lower,
2704 #define FNSZ(upper,lower) &grid_size_ ## lower,
2705
2706 static grid *(*(grid_news[]))(int, int, char*) = { GRIDGEN_LIST(FNNEW) };
2707 static void(*(grid_sizes[]))(int, int, int*, int*, int*) = { GRIDGEN_LIST(FNSZ) };
2708
2709 char *grid_new_desc(grid_type type, int width, int height, random_state *rs)
2710 {
2711     if (type != GRID_PENROSE_P2 && type != GRID_PENROSE_P3)
2712         return NULL;
2713
2714     return grid_new_desc_penrose(type, width, height, rs);
2715 }
2716
2717 char *grid_validate_desc(grid_type type, int width, int height, char *desc)
2718 {
2719     if (type != GRID_PENROSE_P2 && type != GRID_PENROSE_P3) {
2720         if (desc != NULL)
2721             return "Grid description strings not used with this grid type";
2722         return NULL;
2723     }
2724
2725     return grid_validate_desc_penrose(type, width, height, desc);
2726 }
2727
2728 grid *grid_new(grid_type type, int width, int height, char *desc)
2729 {
2730     char *err = grid_validate_desc(type, width, height, desc);
2731     if (err) assert(!"Invalid grid description.");
2732
2733     return grid_news[type](width, height, desc);
2734 }
2735
2736 void grid_compute_size(grid_type type, int width, int height,
2737                        int *tilesize, int *xextent, int *yextent)
2738 {
2739     grid_sizes[type](width, height, tilesize, xextent, yextent);
2740 }
2741
2742 /* ----------- End of grid helpers ------------- */
2743
2744 /* vim: set shiftwidth=4 tabstop=8: */