- 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;
- }
- }
- }
- }
-
- assert(bestdist >= 0);
-
- /* convert to screen coordinates */
- grid_to_screen(ds, g, xbest << shift, ybest << shift,