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 c29da3e56880633bdbcbc67c29e7f25bbedf27d8..20d0ba6c03f401395a27b4eac6ee958e923d6dc9 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.91.2.2 2008/04/11 22:53:16 tgl Exp $
+ *   $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.91.2.3 2009/06/23 16:25:22 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);