Fix an ancient error in dist_ps (distance from point to line segment), which

a number of other geometric operators also depend on.  It miscalculated the
slope of the perpendicular to the given line segment anytime that slope was
other than 0, infinite, or +/-1.  In some cases the error would be masked
because the true closest point on the line segment was one of its endpoints
rather than the intersection point, but in other cases it could give an
arbitrarily bad answer.  Per bug #4872 from Nick Roosevelt.

Bug goes clear back to Berkeley days, so patch all supported branches.
Make a couple of cosmetic adjustments while at it.
This commit is contained in:
Tom Lane 2009-06-23 16:25:02 +00:00
parent 8d355d7bbf
commit b087b018a1

View File

@ -8,7 +8,7 @@
*
*
* IDENTIFICATION
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.101 2009/01/01 17:23:49 momjian Exp $
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.102 2009/06/23 16:25:02 tgl Exp $
*
*-------------------------------------------------------------------------
*/
@ -2449,18 +2449,16 @@ dist_ps_internal(Point *pt, LSEG *lseg)
tmpdist;
Point *ip;
/*
* Construct a line perpendicular to the input segment
* and through the input point
*/
/*
* Construct a line perpendicular to the input segment
* and through the input point
*/
if (lseg->p[1].x == lseg->p[0].x)
m = 0;
else if (lseg->p[1].y == lseg->p[0].y)
{ /* slope is infinite */
m = (double) DBL_MAX;
}
m = (double) DBL_MAX; /* slope is infinite */
else
m = ((lseg->p[0].y - lseg->p[1].y) / (lseg->p[1].x - lseg->p[0].x));
m = (lseg->p[0].x - lseg->p[1].x) / (lseg->p[1].y - lseg->p[0].y);
ln = line_construct_pm(pt, m);
#ifdef GEODEBUG
@ -2469,13 +2467,14 @@ dist_ps_internal(Point *pt, LSEG *lseg)
#endif
/*
* Calculate distance to the line segment or to the endpoints of the
* segment.
* Calculate distance to the line segment or to the nearest endpoint of
* the segment.
*/
/* intersection is on the line segment? */
if ((ip = interpt_sl(lseg, ln)) != NULL)
{
/* yes, so use distance to the intersection point */
result = point_dt(pt, ip);
#ifdef GEODEBUG
printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
@ -2484,7 +2483,7 @@ dist_ps_internal(Point *pt, LSEG *lseg)
}
else
{
/* intersection is not on line segment */
/* no, so use distance to the nearer endpoint */
result = point_dt(pt, &lseg->p[0]);
tmpdist = point_dt(pt, &lseg->p[1]);
if (tmpdist < result)
@ -3790,7 +3789,7 @@ poly_contain(PG_FUNCTION_ARGS)
{
if (point_inside(&(polyb->p[i]), polya->npts, &(polya->p[0])) == 0)
{
#if GEODEBUG
#ifdef GEODEBUG
printf("poly_contain- point (%f,%f) not in polygon\n", polyb->p[i].x, polyb->p[i].y);
#endif
result = false;
@ -3803,7 +3802,7 @@ poly_contain(PG_FUNCTION_ARGS)
{
if (point_inside(&(polya->p[i]), polyb->npts, &(polyb->p[0])) == 1)
{
#if GEODEBUG
#ifdef GEODEBUG
printf("poly_contain- point (%f,%f) in polygon\n", polya->p[i].x, polya->p[i].y);
#endif
result = false;
@ -3814,7 +3813,7 @@ poly_contain(PG_FUNCTION_ARGS)
}
else
{
#if GEODEBUG
#ifdef GEODEBUG
printf("poly_contain- bound box ((%f,%f),(%f,%f)) not inside ((%f,%f),(%f,%f))\n",
polyb->boundbox.low.x, polyb->boundbox.low.y, polyb->boundbox.high.x, polyb->boundbox.high.y,
polya->boundbox.low.x, polya->boundbox.low.y, polya->boundbox.high.x, polya->boundbox.high.y);