+ /*
+ * Now iterate over every point in that bounding box.
+ */
+ xbest = ybest = -1;
+ bestdist = -1;
+ for (y = y0; y <= y1; y++) {
+ for (x = x0; x <= x1; x++) {
+ /*
+ * First, disqualify the point if it's not inside the
+ * polygon, which we work out by counting the edges to the
+ * right of the point. (For tiebreaking purposes when
+ * edges start or end on our y-coordinate or go right
+ * through it, we consider our point to be offset by a
+ * small _positive_ epsilon in both the x- and
+ * y-direction.)
+ */
+ int in = 0;
+ for (i = 0; i < f->order; i++) {
+ int xs = f->edges[i]->dot1->x >> shift;
+ int xe = f->edges[i]->dot2->x >> shift;
+ int ys = f->edges[i]->dot1->y >> shift;
+ int ye = f->edges[i]->dot2->y >> shift;
+ if ((y >= ys && y < ye) || (y >= ye && y < ys)) {
+ /*
+ * The line goes past our y-position. Now we need
+ * to know if its x-coordinate when it does so is
+ * to our right.
+ *
+ * The x-coordinate in question is mathematically
+ * (y - ys) * (xe - xs) / (ye - ys), and we want
+ * to know whether (x - xs) >= that. Of course we
+ * avoid the division, so we can work in integers;
+ * to do this we must multiply both sides of the
+ * inequality by ye - ys, which means we must
+ * first check that's not negative.
+ */
+ int num = xe - xs, denom = ye - ys;
+ if (denom < 0) {
+ num = -num;
+ denom = -denom;
+ }
+ if ((x - xs) * denom >= (y - ys) * num)
+ in ^= 1;
+ }
+ }
+
+ if (in) {
+ long mindist = LONG_MAX;
+
+ /*
+ * This point is inside the polygon, so now we check
+ * its minimum distance to every edge and corner.
+ * First the corners ...
+ */
+ for (i = 0; i < f->order; i++) {
+ int xp = f->dots[i]->x >> shift;
+ int yp = f->dots[i]->y >> shift;
+ int dx = x - xp, dy = y - yp;
+ long dist = (long)dx*dx + (long)dy*dy;
+ if (mindist > dist)
+ mindist = dist;
+ }
+
+ /*
+ * ... and now also check the perpendicular distance
+ * to every edge, if the perpendicular lies between
+ * the edge's endpoints.
+ */
+ for (i = 0; i < f->order; i++) {
+ int xs = f->edges[i]->dot1->x >> shift;
+ int xe = f->edges[i]->dot2->x >> shift;
+ int ys = f->edges[i]->dot1->y >> shift;
+ int ye = f->edges[i]->dot2->y >> shift;
+
+ /*
+ * If s and e are our endpoints, and p our
+ * candidate circle centre, the foot of a
+ * perpendicular from p to the line se lies
+ * between s and e if and only if (p-s).(e-s) lies
+ * strictly between 0 and (e-s).(e-s).
+ */
+ int edx = xe - xs, edy = ye - ys;
+ int pdx = x - xs, pdy = y - ys;
+ long pde = (long)pdx * edx + (long)pdy * edy;
+ long ede = (long)edx * edx + (long)edy * edy;
+ if (0 < pde && pde < ede) {
+ /*
+ * Yes, the nearest point on this edge is
+ * closer than either endpoint, so we must
+ * take it into account by measuring the
+ * perpendicular distance to the edge and
+ * checking its square against mindist.
+ */
+
+ long pdre = (long)pdx * edy - (long)pdy * edx;
+ long sqlen = pdre * pdre / ede;
+
+ if (mindist > sqlen)
+ mindist = sqlen;
+ }
+ }
+
+ /*
+ * Right. Now we know the biggest circle around this
+ * point, so we can check it against bestdist.
+ */
+ if (bestdist < mindist) {
+ bestdist = mindist;
+ xbest = x;
+ ybest = y;
+ }
+ }
+ }