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