mirror of
git://git.tartarus.org/simon/puzzles.git
synced 2025-04-21 16:05:44 -07:00
Replace my brute-force algorithm in face_text_pos with a more complex
but faster and more mathematically sensible one. [originally from svn r9156]
This commit is contained in:
Simon Tatham
566
loopy.c
566
loopy.c
@ -3401,13 +3401,62 @@ static void grid_to_screen(const game_drawstate *ds, const grid *g,
|
|||||||
*y += BORDER(ds->tilesize);
|
*y += BORDER(ds->tilesize);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
static int solve_2x2_matrix(double mx[4], double vin[2], double vout[2])
|
||||||
|
{
|
||||||
|
double inv[4];
|
||||||
|
double det;
|
||||||
|
det = (mx[0]*mx[3] - mx[1]*mx[2]);
|
||||||
|
if (det == 0)
|
||||||
|
return FALSE;
|
||||||
|
|
||||||
|
inv[0] = mx[3] / det;
|
||||||
|
inv[1] = -mx[1] / det;
|
||||||
|
inv[2] = -mx[2] / det;
|
||||||
|
inv[3] = mx[0] / det;
|
||||||
|
|
||||||
|
vout[0] = inv[0]*vin[0] + inv[1]*vin[1];
|
||||||
|
vout[1] = inv[2]*vin[0] + inv[3]*vin[1];
|
||||||
|
|
||||||
|
return TRUE;
|
||||||
|
}
|
||||||
|
|
||||||
|
static int solve_3x3_matrix(double mx[9], double vin[3], double vout[3])
|
||||||
|
{
|
||||||
|
double inv[9];
|
||||||
|
double det;
|
||||||
|
|
||||||
|
det = (mx[0]*mx[4]*mx[8] + mx[1]*mx[5]*mx[6] + mx[2]*mx[3]*mx[7] -
|
||||||
|
mx[0]*mx[5]*mx[7] - mx[1]*mx[3]*mx[8] - mx[2]*mx[4]*mx[6]);
|
||||||
|
if (det == 0)
|
||||||
|
return FALSE;
|
||||||
|
|
||||||
|
inv[0] = (mx[4]*mx[8] - mx[5]*mx[7]) / det;
|
||||||
|
inv[1] = (mx[2]*mx[7] - mx[1]*mx[8]) / det;
|
||||||
|
inv[2] = (mx[1]*mx[5] - mx[2]*mx[4]) / det;
|
||||||
|
inv[3] = (mx[5]*mx[6] - mx[3]*mx[8]) / det;
|
||||||
|
inv[4] = (mx[0]*mx[8] - mx[2]*mx[6]) / det;
|
||||||
|
inv[5] = (mx[2]*mx[3] - mx[0]*mx[5]) / det;
|
||||||
|
inv[6] = (mx[3]*mx[7] - mx[4]*mx[6]) / det;
|
||||||
|
inv[7] = (mx[1]*mx[6] - mx[0]*mx[7]) / det;
|
||||||
|
inv[8] = (mx[0]*mx[4] - mx[1]*mx[3]) / det;
|
||||||
|
|
||||||
|
vout[0] = inv[0]*vin[0] + inv[1]*vin[1] + inv[2]*vin[2];
|
||||||
|
vout[1] = inv[3]*vin[0] + inv[4]*vin[1] + inv[5]*vin[2];
|
||||||
|
vout[2] = inv[6]*vin[0] + inv[7]*vin[1] + inv[8]*vin[2];
|
||||||
|
|
||||||
|
return TRUE;
|
||||||
|
}
|
||||||
|
|
||||||
/* Returns (into x,y) position of centre of face for rendering the text clue.
|
/* Returns (into x,y) position of centre of face for rendering the text clue.
|
||||||
*/
|
*/
|
||||||
static void face_text_pos(const game_drawstate *ds, const grid *g,
|
static void face_text_pos(const game_drawstate *ds, const grid *g,
|
||||||
const grid_face *f, int *xret, int *yret)
|
const grid_face *f, int *xret, int *yret)
|
||||||
{
|
{
|
||||||
int x, y, x0, y0, x1, y1, xbest, ybest, i, shift;
|
double xbest, ybest, bestdist;
|
||||||
long bestdist;
|
int i, j, k, m;
|
||||||
|
grid_dot *edgedot1[3], *edgedot2[3];
|
||||||
|
grid_dot *dots[3];
|
||||||
|
int nedges, ndots;
|
||||||
int faceindex = f - g->faces;
|
int faceindex = f - g->faces;
|
||||||
|
|
||||||
/*
|
/*
|
||||||
@ -3424,151 +3473,422 @@ static void face_text_pos(const game_drawstate *ds, const grid *g,
|
|||||||
* Otherwise, try to find the point in the polygon with the
|
* Otherwise, try to find the point in the polygon with the
|
||||||
* maximum distance to any edge or corner.
|
* maximum distance to any edge or corner.
|
||||||
*
|
*
|
||||||
* Start by working out the face's bounding box, in grid
|
* This point must be in contact with at least three edges and/or
|
||||||
* coordinates.
|
* vertices; so we iterate through all combinations of three of
|
||||||
|
* those, and find candidate points in each set.
|
||||||
|
*
|
||||||
|
* We don't actually iterate literally over _edges_, in the sense
|
||||||
|
* of grid_edge structures. Instead, we fill in edgedot1[] and
|
||||||
|
* edgedot2[] with a pair of dots adjacent in the face's list of
|
||||||
|
* vertices. This ensures that we get the edges in consistent
|
||||||
|
* orientation, which we could not do from the grid structure
|
||||||
|
* alone. (A moment's consideration of an order-3 vertex should
|
||||||
|
* make it clear that if a notional arrow was written on each
|
||||||
|
* edge, _at least one_ of the three faces bordering that vertex
|
||||||
|
* would have to have the two arrows tip-to-tip or tail-to-tail
|
||||||
|
* rather than tip-to-tail.)
|
||||||
*/
|
*/
|
||||||
x0 = x1 = f->dots[0]->x;
|
nedges = ndots = 0;
|
||||||
y0 = y1 = f->dots[0]->y;
|
bestdist = 0;
|
||||||
for (i = 1; i < f->order; i++) {
|
xbest = ybest = 0;
|
||||||
if (x0 > f->dots[i]->x) x0 = f->dots[i]->x;
|
|
||||||
if (x1 < f->dots[i]->x) x1 = f->dots[i]->x;
|
|
||||||
if (y0 > f->dots[i]->y) y0 = f->dots[i]->y;
|
|
||||||
if (y1 < f->dots[i]->y) y1 = f->dots[i]->y;
|
|
||||||
}
|
|
||||||
|
|
||||||
/*
|
for (i = 0; i+2 < 2*f->order; i++) {
|
||||||
* If the grid is at excessive resolution, decide on a scaling
|
if (i < f->order) {
|
||||||
* factor to bring it within reasonable bounds so we don't have to
|
edgedot1[nedges] = f->dots[i];
|
||||||
* think too hard or suffer integer overflow.
|
edgedot2[nedges++] = f->dots[(i+1)%f->order];
|
||||||
*/
|
} else
|
||||||
shift = 0;
|
dots[ndots++] = f->dots[i - f->order];
|
||||||
while (x1 - x0 > 128 || y1 - y0 > 128) {
|
|
||||||
shift++;
|
|
||||||
x0 >>= 1;
|
|
||||||
x1 >>= 1;
|
|
||||||
y0 >>= 1;
|
|
||||||
y1 >>= 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
/*
|
for (j = i+1; j+1 < 2*f->order; j++) {
|
||||||
* Now iterate over every point in that bounding box.
|
if (j < f->order) {
|
||||||
*/
|
edgedot1[nedges] = f->dots[j];
|
||||||
xbest = ybest = -1;
|
edgedot2[nedges++] = f->dots[(j+1)%f->order];
|
||||||
bestdist = -1;
|
} else
|
||||||
for (y = y0; y <= y1; y++) {
|
dots[ndots++] = f->dots[j - f->order];
|
||||||
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) {
|
for (k = j+1; k < 2*f->order; k++) {
|
||||||
long mindist = LONG_MAX;
|
double cx[2], cy[2]; /* candidate positions */
|
||||||
|
int cn = 0; /* number of candidates */
|
||||||
|
|
||||||
|
if (k < f->order) {
|
||||||
|
edgedot1[nedges] = f->dots[k];
|
||||||
|
edgedot2[nedges++] = f->dots[(k+1)%f->order];
|
||||||
|
} else
|
||||||
|
dots[ndots++] = f->dots[k - f->order];
|
||||||
|
|
||||||
/*
|
/*
|
||||||
* This point is inside the polygon, so now we check
|
* Find a point, or pair of points, equidistant from
|
||||||
* its minimum distance to every edge and corner.
|
* all the specified edges and/or vertices.
|
||||||
* First the corners ...
|
|
||||||
*/
|
*/
|
||||||
for (i = 0; i < f->order; i++) {
|
if (nedges == 3) {
|
||||||
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
|
* Three edges. This is a linear matrix equation:
|
||||||
* candidate circle centre, the foot of a
|
* each row of the matrix represents the fact that
|
||||||
* perpendicular from p to the line se lies
|
* the point (x,y) we seek is at distance r from
|
||||||
* between s and e if and only if (p-s).(e-s) lies
|
* that edge, and we solve three of those
|
||||||
* strictly between 0 and (e-s).(e-s).
|
* simultaneously to obtain x,y,r. (We ignore r.)
|
||||||
*/
|
*/
|
||||||
int edx = xe - xs, edy = ye - ys;
|
double matrix[9], vector[3], vector2[3];
|
||||||
int pdx = x - xs, pdy = y - ys;
|
int m;
|
||||||
long pde = (long)pdx * edx + (long)pdy * edy;
|
|
||||||
long ede = (long)edx * edx + (long)edy * edy;
|
for (m = 0; m < 3; m++) {
|
||||||
if (0 < pde && pde < ede) {
|
int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
|
||||||
|
int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
|
||||||
|
int dx = x2-x1, dy = y2-y1;
|
||||||
|
|
||||||
/*
|
/*
|
||||||
* Yes, the nearest point on this edge is
|
* ((x,y) - (x1,y1)) . (dy,-dx) = r |(dx,dy)|
|
||||||
* closer than either endpoint, so we must
|
*
|
||||||
* take it into account by measuring the
|
* => x dy - y dx - r |(dx,dy)| = (x1 dy - y1 dx)
|
||||||
* perpendicular distance to the edge and
|
|
||||||
* checking its square against mindist.
|
|
||||||
*/
|
*/
|
||||||
|
matrix[3*m+0] = dy;
|
||||||
|
matrix[3*m+1] = -dx;
|
||||||
|
matrix[3*m+2] = -sqrt((double)dx*dx+(double)dy*dy);
|
||||||
|
vector[m] = (double)x1*dy - (double)y1*dx;
|
||||||
|
}
|
||||||
|
|
||||||
long pdre = (long)pdx * edy - (long)pdy * edx;
|
if (solve_3x3_matrix(matrix, vector, vector2)) {
|
||||||
long sqlen = pdre * pdre / ede;
|
cx[cn] = vector2[0];
|
||||||
|
cy[cn] = vector2[1];
|
||||||
|
cn++;
|
||||||
|
}
|
||||||
|
} else if (nedges == 2) {
|
||||||
|
/*
|
||||||
|
* Two edges and a dot. This will end up in a
|
||||||
|
* quadratic equation.
|
||||||
|
*
|
||||||
|
* First, look at the two edges. Having our point
|
||||||
|
* be some distance r from both of them gives rise
|
||||||
|
* to a pair of linear equations in x,y,r of the
|
||||||
|
* form
|
||||||
|
*
|
||||||
|
* (x-x1) dy - (y-y1) dx = r sqrt(dx^2+dy^2)
|
||||||
|
*
|
||||||
|
* We eliminate r between those equations to give
|
||||||
|
* us a single linear equation in x,y describing
|
||||||
|
* the locus of points equidistant from both lines
|
||||||
|
* - i.e. the angle bisector.
|
||||||
|
*
|
||||||
|
* We then choose one of x,y to be a parameter t,
|
||||||
|
* and derive linear formulae for x,y,r in terms
|
||||||
|
* of t. This enables us to write down the
|
||||||
|
* circular equation (x-xd)^2+(y-yd)^2=r^2 as a
|
||||||
|
* quadratic in t; solving that and substituting
|
||||||
|
* in for x,y gives us two candidate points.
|
||||||
|
*/
|
||||||
|
double eqs[2][4]; /* a,b,c,d : ax+by+cr=d */
|
||||||
|
double eq[3]; /* a,b,c: ax+by=c */
|
||||||
|
double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
|
||||||
|
double q[3]; /* a,b,c: at^2+bt+c=0 */
|
||||||
|
double disc;
|
||||||
|
|
||||||
if (mindist > sqlen)
|
/* Find equations of the two input lines. */
|
||||||
mindist = sqlen;
|
for (m = 0; m < 2; m++) {
|
||||||
|
int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
|
||||||
|
int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
|
||||||
|
int dx = x2-x1, dy = y2-y1;
|
||||||
|
|
||||||
|
eqs[m][0] = dy;
|
||||||
|
eqs[m][1] = -dx;
|
||||||
|
eqs[m][2] = -sqrt(dx*dx+dy*dy);
|
||||||
|
eqs[m][3] = x1*dy - y1*dx;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Derive the angle bisector by eliminating r. */
|
||||||
|
eq[0] = eqs[0][0]*eqs[1][2] - eqs[1][0]*eqs[0][2];
|
||||||
|
eq[1] = eqs[0][1]*eqs[1][2] - eqs[1][1]*eqs[0][2];
|
||||||
|
eq[2] = eqs[0][3]*eqs[1][2] - eqs[1][3]*eqs[0][2];
|
||||||
|
|
||||||
|
/* Parametrise x and y in terms of some t. */
|
||||||
|
if (abs(eq[0]) < abs(eq[1])) {
|
||||||
|
/* Parameter is x. */
|
||||||
|
xt[0] = 1; xt[1] = 0;
|
||||||
|
yt[0] = -eq[0]/eq[1]; yt[1] = eq[2]/eq[1];
|
||||||
|
} else {
|
||||||
|
/* Parameter is y. */
|
||||||
|
yt[0] = 1; yt[1] = 0;
|
||||||
|
xt[0] = -eq[1]/eq[0]; xt[1] = eq[2]/eq[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Find a linear representation of r using eqs[0]. */
|
||||||
|
rt[0] = -(eqs[0][0]*xt[0] + eqs[0][1]*yt[0])/eqs[0][2];
|
||||||
|
rt[1] = (eqs[0][3] - eqs[0][0]*xt[1] -
|
||||||
|
eqs[0][1]*yt[1])/eqs[0][2];
|
||||||
|
|
||||||
|
/* Construct the quadratic equation. */
|
||||||
|
q[0] = -rt[0]*rt[0];
|
||||||
|
q[1] = -2*rt[0]*rt[1];
|
||||||
|
q[2] = -rt[1]*rt[1];
|
||||||
|
q[0] += xt[0]*xt[0];
|
||||||
|
q[1] += 2*xt[0]*(xt[1]-dots[0]->x);
|
||||||
|
q[2] += (xt[1]-dots[0]->x)*(xt[1]-dots[0]->x);
|
||||||
|
q[0] += yt[0]*yt[0];
|
||||||
|
q[1] += 2*yt[0]*(yt[1]-dots[0]->y);
|
||||||
|
q[2] += (yt[1]-dots[0]->y)*(yt[1]-dots[0]->y);
|
||||||
|
|
||||||
|
/* And solve it. */
|
||||||
|
disc = q[1]*q[1] - 4*q[0]*q[2];
|
||||||
|
if (disc >= 0) {
|
||||||
|
double t;
|
||||||
|
|
||||||
|
disc = sqrt(disc);
|
||||||
|
|
||||||
|
t = (-q[1] + disc) / (2*q[0]);
|
||||||
|
cx[cn] = xt[0]*t + xt[1];
|
||||||
|
cy[cn] = yt[0]*t + yt[1];
|
||||||
|
cn++;
|
||||||
|
|
||||||
|
t = (-q[1] - disc) / (2*q[0]);
|
||||||
|
cx[cn] = xt[0]*t + xt[1];
|
||||||
|
cy[cn] = yt[0]*t + yt[1];
|
||||||
|
cn++;
|
||||||
|
}
|
||||||
|
} else if (nedges == 1) {
|
||||||
|
/*
|
||||||
|
* Two dots and an edge. This one's another
|
||||||
|
* quadratic equation.
|
||||||
|
*
|
||||||
|
* The point we want must lie on the perpendicular
|
||||||
|
* bisector of the two dots; that much is obvious.
|
||||||
|
* So we can construct a parametrisation of that
|
||||||
|
* bisecting line, giving linear formulae for x,y
|
||||||
|
* in terms of t. We can also express the distance
|
||||||
|
* from the edge as such a linear formula.
|
||||||
|
*
|
||||||
|
* Then we set that equal to the radius of the
|
||||||
|
* circle passing through the two points, which is
|
||||||
|
* a Pythagoras exercise; that gives rise to a
|
||||||
|
* quadratic in t, which we solve.
|
||||||
|
*/
|
||||||
|
double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
|
||||||
|
double q[3]; /* a,b,c: at^2+bt+c=0 */
|
||||||
|
double disc;
|
||||||
|
double halfsep;
|
||||||
|
|
||||||
|
/* Find parametric formulae for x,y. */
|
||||||
|
{
|
||||||
|
int x1 = dots[0]->x, x2 = dots[1]->x;
|
||||||
|
int y1 = dots[0]->y, y2 = dots[1]->y;
|
||||||
|
int dx = x2-x1, dy = y2-y1;
|
||||||
|
double d = sqrt((double)dx*dx + (double)dy*dy);
|
||||||
|
|
||||||
|
xt[1] = (x1+x2)/2.0;
|
||||||
|
yt[1] = (y1+y2)/2.0;
|
||||||
|
/* It's convenient if we have t at standard scale. */
|
||||||
|
xt[0] = -dy/d;
|
||||||
|
yt[0] = dx/d;
|
||||||
|
|
||||||
|
/* Also note down half the separation between
|
||||||
|
* the dots, for use in computing the circle radius. */
|
||||||
|
halfsep = 0.5*d;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Find a parametric formula for r. */
|
||||||
|
{
|
||||||
|
int x1 = edgedot1[0]->x, x2 = edgedot2[0]->x;
|
||||||
|
int y1 = edgedot1[0]->y, y2 = edgedot2[0]->y;
|
||||||
|
int dx = x2-x1, dy = y2-y1;
|
||||||
|
double d = sqrt((double)dx*dx + (double)dy*dy);
|
||||||
|
rt[0] = (xt[0]*dy - yt[0]*dx) / d;
|
||||||
|
rt[1] = ((xt[1]-x1)*dy - (yt[1]-y1)*dx) / d;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Construct the quadratic equation. */
|
||||||
|
q[0] = rt[0]*rt[0];
|
||||||
|
q[1] = 2*rt[0]*rt[1];
|
||||||
|
q[2] = rt[1]*rt[1];
|
||||||
|
q[0] -= 1;
|
||||||
|
q[2] -= halfsep*halfsep;
|
||||||
|
|
||||||
|
/* And solve it. */
|
||||||
|
disc = q[1]*q[1] - 4*q[0]*q[2];
|
||||||
|
if (disc >= 0) {
|
||||||
|
double t;
|
||||||
|
|
||||||
|
disc = sqrt(disc);
|
||||||
|
|
||||||
|
t = (-q[1] + disc) / (2*q[0]);
|
||||||
|
cx[cn] = xt[0]*t + xt[1];
|
||||||
|
cy[cn] = yt[0]*t + yt[1];
|
||||||
|
cn++;
|
||||||
|
|
||||||
|
t = (-q[1] - disc) / (2*q[0]);
|
||||||
|
cx[cn] = xt[0]*t + xt[1];
|
||||||
|
cy[cn] = yt[0]*t + yt[1];
|
||||||
|
cn++;
|
||||||
|
}
|
||||||
|
} else if (nedges == 0) {
|
||||||
|
/*
|
||||||
|
* Three dots. This is another linear matrix
|
||||||
|
* equation, this time with each row of the matrix
|
||||||
|
* representing the perpendicular bisector between
|
||||||
|
* two of the points. Of course we only need two
|
||||||
|
* such lines to find their intersection, so we
|
||||||
|
* need only solve a 2x2 matrix equation.
|
||||||
|
*/
|
||||||
|
|
||||||
|
double matrix[4], vector[2], vector2[2];
|
||||||
|
int m;
|
||||||
|
|
||||||
|
for (m = 0; m < 2; m++) {
|
||||||
|
int x1 = dots[m]->x, x2 = dots[m+1]->x;
|
||||||
|
int y1 = dots[m]->y, y2 = dots[m+1]->y;
|
||||||
|
int dx = x2-x1, dy = y2-y1;
|
||||||
|
|
||||||
|
/*
|
||||||
|
* ((x,y) - (x1,y1)) . (dx,dy) = 1/2 |(dx,dy)|^2
|
||||||
|
*
|
||||||
|
* => 2x dx + 2y dy = dx^2+dy^2 + (2 x1 dx + 2 y1 dy)
|
||||||
|
*/
|
||||||
|
matrix[2*m+0] = 2*dx;
|
||||||
|
matrix[2*m+1] = 2*dy;
|
||||||
|
vector[m] = ((double)dx*dx + (double)dy*dy +
|
||||||
|
2.0*x1*dx + 2.0*y1*dy);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (solve_2x2_matrix(matrix, vector, vector2)) {
|
||||||
|
cx[cn] = vector2[0];
|
||||||
|
cy[cn] = vector2[1];
|
||||||
|
cn++;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/*
|
/*
|
||||||
* Right. Now we know the biggest circle around this
|
* Now go through our candidate points and see if any
|
||||||
* point, so we can check it against bestdist.
|
* of them are better than what we've got so far.
|
||||||
*/
|
*/
|
||||||
if (bestdist < mindist) {
|
for (m = 0; m < cn; m++) {
|
||||||
bestdist = mindist;
|
double x = cx[m], y = cy[m];
|
||||||
xbest = x;
|
|
||||||
ybest = y;
|
/*
|
||||||
|
* 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 e, in = 0;
|
||||||
|
for (e = 0; e < f->order; e++) {
|
||||||
|
int xs = f->edges[e]->dot1->x;
|
||||||
|
int xe = f->edges[e]->dot2->x;
|
||||||
|
int ys = f->edges[e]->dot1->y;
|
||||||
|
int ye = f->edges[e]->dot2->y;
|
||||||
|
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) {
|
||||||
|
double mindist = HUGE_VAL;
|
||||||
|
int e, d;
|
||||||
|
|
||||||
|
/*
|
||||||
|
* This point is inside the polygon, so now we check
|
||||||
|
* its minimum distance to every edge and corner.
|
||||||
|
* First the corners ...
|
||||||
|
*/
|
||||||
|
for (d = 0; d < f->order; d++) {
|
||||||
|
int xp = f->dots[d]->x;
|
||||||
|
int yp = f->dots[d]->y;
|
||||||
|
double dx = x - xp, dy = y - yp;
|
||||||
|
double dist = dx*dx + 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 (e = 0; e < f->order; e++) {
|
||||||
|
int xs = f->edges[e]->dot1->x;
|
||||||
|
int xe = f->edges[e]->dot2->x;
|
||||||
|
int ys = f->edges[e]->dot1->y;
|
||||||
|
int ye = f->edges[e]->dot2->y;
|
||||||
|
|
||||||
|
/*
|
||||||
|
* 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;
|
||||||
|
double pdx = x - xs, pdy = y - ys;
|
||||||
|
double pde = pdx * edx + 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.
|
||||||
|
*/
|
||||||
|
|
||||||
|
double pdre = pdx * edy - pdy * edx;
|
||||||
|
double 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;
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
if (k < f->order)
|
||||||
|
nedges--;
|
||||||
|
else
|
||||||
|
ndots--;
|
||||||
}
|
}
|
||||||
|
if (j < f->order)
|
||||||
|
nedges--;
|
||||||
|
else
|
||||||
|
ndots--;
|
||||||
}
|
}
|
||||||
|
if (i < f->order)
|
||||||
|
nedges--;
|
||||||
|
else
|
||||||
|
ndots--;
|
||||||
}
|
}
|
||||||
|
|
||||||
assert(bestdist >= 0);
|
assert(bestdist > 0);
|
||||||
|
|
||||||
/* convert to screen coordinates */
|
/* convert to screen coordinates. Round doubles to nearest. */
|
||||||
grid_to_screen(ds, g, xbest << shift, ybest << shift,
|
grid_to_screen(ds, g, xbest+0.5, ybest+0.5,
|
||||||
&ds->textx[faceindex], &ds->texty[faceindex]);
|
&ds->textx[faceindex], &ds->texty[faceindex]);
|
||||||
|
|
||||||
*xret = ds->textx[faceindex];
|
*xret = ds->textx[faceindex];
|
||||||
|
|||||||
Reference in New Issue
Block a user