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