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.

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index c2444c38e889b1227f095c4d3ba49c2e8427cdfc..24157a53c59955f70abb2669df8f3fb57707d3da 100644 (file)
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -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);