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