Fix problems in handling the line data type

According to the source history, the internal format of line data type
has changed, but various functions working with it did were not updated
and thus were producing wrong results.

This patch addresses various such issues, in particular:

* Reject invalid specification A=B=0 on receive
* Reject same points on line_construct_pp()
* Fix perpendicular operator when negative values are involved
* Avoid division by zero on perpendicular operator
* Fix intersection and distance operators when neither A nor B are 1
* Return NULL for closest point when objects are parallel
* Check whether closest point of line segments is the intersection point
* Fix closest point of line segments being on the wrong segment

Aside from handling those issues, the patch also aims to make operators
more symmetric and less sen to precision loss.  The EPSILON interferes
with even minor changes, but the least we can do is applying it to both
sides of the operators equally.

Author: Emre Hasegeli
Reviewed-by: Tomas Vondra
Discussion: https://www.postgresql.org/message-id/CAE2gYzxF7-5djV6-cEvqQu-fNsnt%3DEqbOURx7ZDg%2BVv6ZMTWbg%40mail.gmail.com

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 13dc1ab6e6002e505df49560b19fe3bda52a4c38..e7c1160131fa70b8446fa0cdf75b1444f6e28907 100644 (file)
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -48,6 +48,7 @@ static int    point_inside(Point *p, int npts, Point *plist);
 
 /* Routines for lines */
 static inline void line_construct(LINE *result, Point *pt, float8 m);
+static inline float8 line_sl(LINE *line);
 static inline float8 line_invsl(LINE *line);
 static bool    line_interpt_line(Point *result, LINE *l1, LINE *l2);
 static bool line_contain_point(LINE *line, Point *point);
@@ -980,6 +981,11 @@ line_recv(PG_FUNCTION_ARGS)
        line->B = pq_getmsgfloat8(buf);
        line->C = pq_getmsgfloat8(buf);
 
+       if (FPzero(line->A) && FPzero(line->B))
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
+                                errmsg("invalid line specification: A and B cannot both be zero")));
+
        PG_RETURN_LINE_P(line);
 }
 
@@ -1040,6 +1046,11 @@ line_construct_pp(PG_FUNCTION_ARGS)
        Point      *pt2 = PG_GETARG_POINT_P(1);
        LINE       *result = (LINE *) palloc(sizeof(LINE));
 
+       if (point_eq_point(pt1, pt2))
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+                                errmsg("invalid line specification: must be two distinct points")));
+
        line_construct(result, pt1, point_sl(pt1, pt2));
 
        PG_RETURN_LINE_P(result);
@@ -1076,11 +1087,15 @@ line_perp(PG_FUNCTION_ARGS)
 
        if (FPzero(l1->A))
                PG_RETURN_BOOL(FPzero(l2->B));
-       else if (FPzero(l1->B))
+       if (FPzero(l2->A))
+               PG_RETURN_BOOL(FPzero(l1->B));
+       if (FPzero(l1->B))
                PG_RETURN_BOOL(FPzero(l2->A));
+       if (FPzero(l2->B))
+               PG_RETURN_BOOL(FPzero(l1->A));
 
-       PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B),
-                                                                  float8_mul(l1->B, l2->A)), -1.0));
+       PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->A),
+                                                                  float8_mul(l1->B, l2->B)), -1.0));
 }
 
 Datum
@@ -1135,6 +1150,20 @@ line_eq(PG_FUNCTION_ARGS)
  *     Line arithmetic routines.
  *---------------------------------------------------------*/
 
+/*
+ * Return slope of the line
+ */
+static inline float8
+line_sl(LINE *line)
+{
+       if (FPzero(line->A))
+               return 0.0;
+       if (FPzero(line->B))
+               return DBL_MAX;
+       return float8_div(line->A, -line->B);
+}
+
+
 /*
  * Return inverse slope of the line
  */
@@ -1157,16 +1186,21 @@ line_distance(PG_FUNCTION_ARGS)
 {
        LINE       *l1 = PG_GETARG_LINE_P(0);
        LINE       *l2 = PG_GETARG_LINE_P(1);
-       float8          result;
-       Point           tmp;
+       float8          ratio;
 
        if (line_interpt_line(NULL, l1, l2))    /* intersecting? */
                PG_RETURN_FLOAT8(0.0);
-       if (FPzero(l1->B))                      /* vertical? */
-               PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C)));
-       point_construct(&tmp, 0.0, l1->C);
-       result = line_closept_point(NULL, l2, &tmp);
-       PG_RETURN_FLOAT8(result);
+
+       if (!FPzero(l1->A) && !isnan(l1->A) && !FPzero(l2->A) && !isnan(l2->A))
+               ratio = float8_div(l1->A, l2->A);
+       else if (!FPzero(l1->B) && !isnan(l1->B) && !FPzero(l2->B) && !isnan(l2->B))
+               ratio = float8_div(l1->B, l2->B);
+       else
+               ratio = 1.0;
+
+       PG_RETURN_FLOAT8(float8_div(fabs(float8_mi(l1->C,
+                                                                                          float8_mul(ratio, l2->C))),
+                                                               HYPOT(l1->A, l1->B)));
 }
 
 /* line_interpt()
@@ -1207,27 +1241,36 @@ line_interpt_line(Point *result, LINE *l1, LINE *l2)
        float8          x,
                                y;
 
-       if (FPzero(l1->B))                      /* l1 vertical? */
+       if (!FPzero(l1->B))
        {
-               if (FPzero(l2->B))              /* l2 vertical? */
+               if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
                        return false;
 
-               x = l1->C;
-               y = float8_pl(float8_mul(l2->A, x), l2->C);
-       }
-       else if (FPzero(l2->B))         /* l2 vertical? */
-       {
-               x = l2->C;
-               y = float8_pl(float8_mul(l1->A, x), l1->C);
+               x = float8_div(float8_mi(float8_mul(l1->B, l2->C),
+                                                                float8_mul(l2->B, l1->C)),
+                                          float8_mi(float8_mul(l1->A, l2->B),
+                                                                float8_mul(l2->A, l1->B)));
+               y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B);
        }
-       else
+       else if (!FPzero(l2->B))
        {
-               if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
+               if (FPeq(l1->A, float8_mul(l2->A, float8_div(l1->B, l2->B))))
                        return false;
 
-               x = float8_div(float8_mi(l1->C, l2->C), float8_mi(l2->A, l1->A));
-               y = float8_pl(float8_mul(l1->A, x), l1->C);
+               x = float8_div(float8_mi(float8_mul(l2->B, l1->C),
+                                                                float8_mul(l1->B, l2->C)),
+                                          float8_mi(float8_mul(l2->A, l1->B),
+                                                                float8_mul(l1->A, l2->B)));
+               y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B);
        }
+       else
+               return false;
+
+       /* On some platforms, the preceding expressions tend to produce -0. */
+       if (x == 0.0)
+               x = 0.0;
+       if (y == 0.0)
+               y = 0.0;
 
        if (result != NULL)
                point_construct(result, x, y);
@@ -2347,16 +2390,8 @@ dist_sl(PG_FUNCTION_ARGS)
 {
        LSEG       *lseg = PG_GETARG_LSEG_P(0);
        LINE       *line = PG_GETARG_LINE_P(1);
-       float8          result;
-
-       if (lseg_interpt_line(NULL, lseg, line))
-               result = 0.0;
-       else
-               /* XXX shouldn't we take the min not max? */
-               result = float8_max(line_closept_point(NULL, line, &lseg->p[0]),
-                                                       line_closept_point(NULL, line, &lseg->p[1]));
 
-       PG_RETURN_FLOAT8(result);
+       PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line));
 }
 
 /*
@@ -2533,20 +2568,26 @@ lseg_interpt_line(Point *result, LSEG *lseg, LINE *line)
 static float8
 line_closept_point(Point *result, LINE *line, Point *point)
 {
-       bool            retval PG_USED_FOR_ASSERTS_ONLY;
-       Point       closept;
+       Point           closept;
        LINE            tmp;
 
-       /* We drop a perpendicular to find the intersection point. */
+       /*
+        * We drop a perpendicular to find the intersection point.  Ordinarily
+        * we should always find it, but that can fail in the presence of NaN
+        * coordinates, and perhaps even from simple roundoff issues.
+        */
        line_construct(&tmp, point, line_invsl(line));
-       retval = line_interpt_line(&closept, line, &tmp);
+       if (!line_interpt_line(&closept, &tmp, line))
+       {
+               if (result != NULL)
+                       *result = *point;
 
-       Assert(retval); /* perpendicular lines always intersect */
+               return get_float8_nan();
+       }
 
        if (result != NULL)
                *result = closept;
 
-       /* Then we calculate the distance between the points. */
        return point_dt(&closept, point);
 }
 
@@ -2620,27 +2661,38 @@ lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2)
        float8          dist,
                                d;
 
-       d = lseg_closept_point(NULL, l1, &l2->p[0]);
-       dist = d;
-       if (result != NULL)
-               *result = l2->p[0];
+       /* First, we handle the case when the line segments are intersecting. */
+       if (lseg_interpt_lseg(result, l1, l2))
+               return 0.0;
 
-       d = lseg_closept_point(NULL, l1, &l2->p[1]);
+       /*
+        * Then, we find the closest points from the endpoints of the second
+        * line segment, and keep the closest one.
+        */
+       dist = lseg_closept_point(result, l1, &l2->p[0]);
+       d = lseg_closept_point(&point, l1, &l2->p[1]);
        if (float8_lt(d, dist))
        {
                dist = d;
                if (result != NULL)
-                       *result = l2->p[1];
+                       *result = point;
        }
 
-       if (float8_lt(lseg_closept_point(&point, l2, &l1->p[0]), dist))
-               d = lseg_closept_point(result, l1, &point);
-
-       if (float8_lt(lseg_closept_point(&point, l2, &l1->p[1]), dist))
-               d = lseg_closept_point(result, l1, &point);
-
+       /* The closest point can still be one of the endpoints, so we test them. */
+       d = lseg_closept_point(NULL, l2, &l1->p[0]);
+       if (float8_lt(d, dist))
+       {
+               dist = d;
+               if (result != NULL)
+                       *result = l1->p[0];
+       }
+       d = lseg_closept_point(NULL, l2, &l1->p[1]);
        if (float8_lt(d, dist))
+       {
                dist = d;
+               if (result != NULL)
+                       *result = l1->p[1];
+       }
 
        return dist;
 }
@@ -2652,6 +2704,9 @@ close_lseg(PG_FUNCTION_ARGS)
        LSEG       *l2 = PG_GETARG_LSEG_P(1);
        Point      *result;
 
+       if (lseg_sl(l1) == lseg_sl(l2))
+               PG_RETURN_NULL();
+
        result = (Point *) palloc(sizeof(Point));
 
        if (isnan(lseg_closept_lseg(result, l2, l1)))
@@ -2826,6 +2881,9 @@ close_ls(PG_FUNCTION_ARGS)
        LSEG       *lseg = PG_GETARG_LSEG_P(1);
        Point      *result;
 
+       if (lseg_sl(lseg) == line_sl(line))
+               PG_RETURN_NULL();
+
        result = (Point *) palloc(sizeof(Point));
 
        if (isnan(lseg_closept_line(result, lseg, line)))