+ /*
+ * Analyse the map to find a canonical line segment
+ * corresponding to each edge, and a canonical point
+ * corresponding to each region. The former are where we'll
+ * eventually put error markers; the latter are where we'll put
+ * per-region flags such as numbers (when in diagnostic mode).
+ */
+ {
+ int *bestx, *besty, *an, pass;
+ float *ax, *ay, *best;
+
+ ax = snewn(state->map->ngraph + n, float);
+ ay = snewn(state->map->ngraph + n, float);
+ an = snewn(state->map->ngraph + n, int);
+ bestx = snewn(state->map->ngraph + n, int);
+ besty = snewn(state->map->ngraph + n, int);
+ best = snewn(state->map->ngraph + n, float);
+
+ for (i = 0; i < state->map->ngraph + n; i++) {
+ bestx[i] = besty[i] = -1;
+ best[i] = (float)(2*(w+h)+1);
+ ax[i] = ay[i] = 0.0F;
+ an[i] = 0;
+ }
+
+ /*
+ * We make two passes over the map, finding all the line
+ * segments separating regions and all the suitable points
+ * within regions. In the first pass, we compute the
+ * _average_ x and y coordinate of all the points in a
+ * given class; in the second pass, for each such average
+ * point, we find the candidate closest to it and call that
+ * canonical.
+ *
+ * Line segments are considered to have coordinates in
+ * their centre. Thus, at least one coordinate for any line
+ * segment is always something-and-a-half; so we store our
+ * coordinates as twice their normal value.
+ */
+ for (pass = 0; pass < 2; pass++) {
+ int x, y;
+
+ for (y = 0; y < h; y++)
+ for (x = 0; x < w; x++) {
+ int ex[4], ey[4], ea[4], eb[4], en = 0;
+
+ /*
+ * Look for an edge to the right of this
+ * square, an edge below it, and an edge in the
+ * middle of it. Also look to see if the point
+ * at the bottom right of this square is on an
+ * edge (and isn't a place where more than two
+ * regions meet).
+ */
+ if (x+1 < w) {
+ /* right edge */
+ ea[en] = state->map->map[RE * wh + y*w+x];
+ eb[en] = state->map->map[LE * wh + y*w+(x+1)];
+ ex[en] = (x+1)*2;
+ ey[en] = y*2+1;
+ en++;
+ }
+ if (y+1 < h) {
+ /* bottom edge */
+ ea[en] = state->map->map[BE * wh + y*w+x];
+ eb[en] = state->map->map[TE * wh + (y+1)*w+x];
+ ex[en] = x*2+1;
+ ey[en] = (y+1)*2;
+ en++;
+ }
+ /* diagonal edge */
+ ea[en] = state->map->map[TE * wh + y*w+x];
+ eb[en] = state->map->map[BE * wh + y*w+x];
+ ex[en] = x*2+1;
+ ey[en] = y*2+1;
+ en++;
+
+ if (x+1 < w && y+1 < h) {
+ /* bottom right corner */
+ int oct[8], othercol, nchanges;
+ oct[0] = state->map->map[RE * wh + y*w+x];
+ oct[1] = state->map->map[LE * wh + y*w+(x+1)];
+ oct[2] = state->map->map[BE * wh + y*w+(x+1)];
+ oct[3] = state->map->map[TE * wh + (y+1)*w+(x+1)];
+ oct[4] = state->map->map[LE * wh + (y+1)*w+(x+1)];
+ oct[5] = state->map->map[RE * wh + (y+1)*w+x];
+ oct[6] = state->map->map[TE * wh + (y+1)*w+x];
+ oct[7] = state->map->map[BE * wh + y*w+x];
+
+ othercol = -1;
+ nchanges = 0;
+ for (i = 0; i < 8; i++) {
+ if (oct[i] != oct[0]) {
+ if (othercol < 0)
+ othercol = oct[i];
+ else if (othercol != oct[i])
+ break; /* three colours at this point */
+ }
+ if (oct[i] != oct[(i+1) & 7])
+ nchanges++;
+ }
+
+ /*
+ * Now if there are exactly two regions at
+ * this point (not one, and not three or
+ * more), and only two changes around the
+ * loop, then this is a valid place to put
+ * an error marker.
+ */
+ if (i == 8 && othercol >= 0 && nchanges == 2) {
+ ea[en] = oct[0];
+ eb[en] = othercol;
+ ex[en] = (x+1)*2;
+ ey[en] = (y+1)*2;
+ en++;
+ }
+
+ /*
+ * If there's exactly _one_ region at this
+ * point, on the other hand, it's a valid
+ * place to put a region centre.
+ */
+ if (othercol < 0) {
+ ea[en] = eb[en] = oct[0];
+ ex[en] = (x+1)*2;
+ ey[en] = (y+1)*2;
+ en++;
+ }
+ }
+
+ /*
+ * Now process the points we've found, one by
+ * one.
+ */
+ for (i = 0; i < en; i++) {
+ int emin = min(ea[i], eb[i]);
+ int emax = max(ea[i], eb[i]);
+ int gindex;
+
+ if (emin != emax) {
+ /* Graph edge */
+ gindex =
+ graph_edge_index(state->map->graph, n,
+ state->map->ngraph, emin,
+ emax);
+ } else {
+ /* Region number */
+ gindex = state->map->ngraph + emin;
+ }
+
+ assert(gindex >= 0);
+
+ if (pass == 0) {
+ /*
+ * In pass 0, accumulate the values
+ * we'll use to compute the average
+ * positions.
+ */
+ ax[gindex] += ex[i];
+ ay[gindex] += ey[i];
+ an[gindex] += 1;
+ } else {
+ /*
+ * In pass 1, work out whether this
+ * point is closer to the average than
+ * the last one we've seen.
+ */
+ float dx, dy, d;
+
+ assert(an[gindex] > 0);
+ dx = ex[i] - ax[gindex];
+ dy = ey[i] - ay[gindex];
+ d = (float)sqrt(dx*dx + dy*dy);
+ if (d < best[gindex]) {
+ best[gindex] = d;
+ bestx[gindex] = ex[i];
+ besty[gindex] = ey[i];
+ }
+ }
+ }
+ }
+
+ if (pass == 0) {
+ for (i = 0; i < state->map->ngraph + n; i++)
+ if (an[i] > 0) {
+ ax[i] /= an[i];
+ ay[i] /= an[i];
+ }
+ }
+ }
+
+ state->map->edgex = snewn(state->map->ngraph, int);
+ state->map->edgey = snewn(state->map->ngraph, int);
+ memcpy(state->map->edgex, bestx, state->map->ngraph * sizeof(int));
+ memcpy(state->map->edgey, besty, state->map->ngraph * sizeof(int));
+
+ state->map->regionx = snewn(n, int);
+ state->map->regiony = snewn(n, int);
+ memcpy(state->map->regionx, bestx + state->map->ngraph, n*sizeof(int));
+ memcpy(state->map->regiony, besty + state->map->ngraph, n*sizeof(int));
+
+ for (i = 0; i < state->map->ngraph; i++)
+ if (state->map->edgex[i] < 0) {
+ /* Find the other representation of this edge. */
+ int e = state->map->graph[i];
+ int iprime = graph_edge_index(state->map->graph, n,
+ state->map->ngraph, e%n, e/n);
+ assert(state->map->edgex[iprime] >= 0);
+ state->map->edgex[i] = state->map->edgex[iprime];
+ state->map->edgey[i] = state->map->edgey[iprime];
+ }
+
+ sfree(ax);
+ sfree(ay);
+ sfree(an);
+ sfree(best);
+ sfree(bestx);
+ sfree(besty);
+ }
+