From 82fd7f782555f7390210e8c07f811bb1ca419659 Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sun, 28 May 2017 11:35:17 +0200 Subject: [PATCH 4/4] line-fixes-v03 Fix obvious problems around the line datatype I have noticed some line operators retuning wrong results, and Tom Lane spotted similar problems on more places. Source history reveals that during 1990s, the internal format of the line datatype is changed, but most functions haven't got the hint. The fixes include: * Reject invalid specification A=B=0 on receive * Avoid division by zero on perpendicular operator * Fix intersection and distance operators when neither A nor B are 1 * Avoid point distance operator crashing on float precision loss * Fix line segment distance by getting the minimum as the comment suggested The changes are also aiming make line operators more symmetric. Changing results for same lines in different order are surely not expected. Previous discussion: https://www.postgresql.org/message-id/flat/CAE2gYzw_-z%3DV2kh8QqFjenu%3D8MJXzOP44wRW%3DAzzeamrmTT1%3DQ%40mail.gmail.com --- src/backend/utils/adt/geo_ops.c | 134 ++++++++++++++++++++++++---------------- 1 file changed, 81 insertions(+), 53 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index dbf1de0ca2..c8714ff78f 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -1003,20 +1003,25 @@ line_recv(PG_FUNCTION_ARGS) { StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); LINE *line; line = (LINE *) palloc(sizeof(LINE)); line->A = pq_getmsgfloat8(buf); line->B = pq_getmsgfloat8(buf); line->C = pq_getmsgfloat8(buf); + if (line->A == 0.0 && line->B == 0.0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid line specification: A and B cannot both be zero"))); + PG_RETURN_LINE_P(line); } /* * line_send - converts line to binary format */ Datum line_send(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); @@ -1123,25 +1128,28 @@ line_parallel(PG_FUNCTION_ARGS) } Datum line_perp(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); if (FPzero(l1->A)) PG_RETURN_BOOL(FPzero(l2->B)); - else if (FPzero(l1->B)) + if (FPzero(l1->B)) PG_RETURN_BOOL(FPzero(l2->A)); + if (FPzero(l2->A)) + PG_RETURN_BOOL(FPzero(l1->B)); + 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_mul(l1->A, l2->B), -float8_mul(l1->B, l2->A))); } Datum line_vertical(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); PG_RETURN_BOOL(FPzero(line->B)); } @@ -1153,20 +1161,21 @@ line_horizontal(PG_FUNCTION_ARGS) PG_RETURN_BOOL(FPzero(line->A)); } Datum line_eq(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); float8 ratio; + /* A and B cannot be both 0, but we never really know with the EPSILON. */ if (!FPzero(l2->A)) ratio = float8_div(l1->A, l2->A); else if (!FPzero(l2->B)) ratio = float8_div(l1->B, l2->B); else if (!FPzero(l2->C)) ratio = float8_div(l1->C, l2->C); else ratio = 1.0; PG_RETURN_BOOL(FPeq(l1->A, float8_mul(ratio, l2->A)) && @@ -1180,30 +1189,36 @@ line_eq(PG_FUNCTION_ARGS) *---------------------------------------------------------*/ /* line_distance() * Distance between two lines. */ Datum 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_internal(NULL, l1, l2)) 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 = dist_pl_internal(&tmp, l2); - PG_RETURN_FLOAT8(result); + + /* A and B cannot be both 0, but we never really know with the EPSILON. */ + if (!FPzero(l2->A)) + ratio = float8_div(l1->A, l2->A); + else if (!FPzero(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))), + float8_hypot(l1->A, l1->B))); } /* line_interpt() * Point where two lines l1, l2 intersect (if any) */ Datum line_interpt(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); @@ -1226,38 +1241,62 @@ line_interpt(PG_FUNCTION_ARGS) * NOTE: If the lines are identical then we will find they are parallel * and report "no intersection". This is a little weird, but since * there's no *unique* intersection, maybe it's appropriate behavior. */ static bool line_interpt_internal(Point *result, LINE *l1, LINE *l2) { float8 x, y; - if (FPzero(l1->B)) /* l1 vertical? */ + if (FPzero(l1->A)) /* l1 horizontal? */ + { + if (FPzero(l2->A)) /* l2 horizontal? */ + return false; + + y = float8_div(-l1->C, l1->B); + x = float8_div(-float8_pl(float8_mul(l2->B, y), l2->C), l2->A); + } + else if (FPzero(l1->B)) /* l1 vertical? */ { if (FPzero(l2->B)) /* l2 vertical? */ return false; - x = l1->C; - y = float8_pl(float8_mul(l2->A, x), l2->C); + x = float8_div(-l1->C, l1->A); + y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B); + } + else if (FPzero(l2->A)) /* l2 horizontal? */ + { + if (FPzero(l1->A)) /* l1 horizontal? */ + return false; + + y = float8_div(-l2->C, l2->B); + x = float8_div(-float8_pl(float8_mul(l1->B, y), l1->C), l1->A); + } + else if (FPzero(l2->B)) /* l2 vertical? */ + { + if (FPzero(l1->B)) /* l1 vertical? */ + return false; + + x = float8_div(-l2->C, l2->A); + y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B); } else { - if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B)))) + if (FPeq(float8_mul(l1->A, l2->B), float8_mul(l2->A, l1->B))) return false; - if (FPzero(l2->B)) - x = l2->C; - else - 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(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); } if (result != NULL) point_construct(result, x, y); #ifdef GEODEBUG printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n", DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C); printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y); #endif @@ -2061,46 +2100,35 @@ lseg_parallel(PG_FUNCTION_ARGS) m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y); PG_RETURN_BOOL(FPeq(m1, m2)); } /* lseg_perp() * Determine if two line segments are perpendicular. * - * This code did not get the correct answer for - * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg - * So, modified it to check explicitly for slope of vertical line - * returned by float8_slope() and the results seem better. - * - thomas 1998-01-31 + * We are comparing the slope of the first line segment with the inverse slope + * of the second line segment. */ Datum lseg_perp(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); float8 m1, m2; m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); - m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y); + m2 = float8_slope(l2->p[0].y, l2->p[1].y, l2->p[1].x, l2->p[0].x); -#ifdef GEODEBUG - printf("lseg_perp- slopes are %g and %g\n", m1, m2); -#endif - if (FPzero(m1)) - PG_RETURN_BOOL(FPeq(m2, DBL_MAX)); - else if (FPzero(m2)) - PG_RETURN_BOOL(FPeq(m1, DBL_MAX)); - - PG_RETURN_BOOL(FPeq(float8_div(m1, m2), -1.0)); + PG_RETURN_BOOL(FPeq(m1, m2)); } Datum lseg_vertical(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x)); } @@ -2467,22 +2495,21 @@ dist_pb(PG_FUNCTION_ARGS) Datum 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_internal(NULL, lseg, line)) result = 0.0; else - /* XXX shouldn't we take the min not max? */ - result = float8_max(dist_pl_internal(&lseg->p[0], line), + result = float8_min(dist_pl_internal(&lseg->p[0], line), dist_pl_internal(&lseg->p[1], line)); PG_RETURN_FLOAT8(result); } /* * Distance from a lseg to a box */ Datum dist_sb(PG_FUNCTION_ARGS) @@ -2663,43 +2690,44 @@ lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line) /* close_pl - * The intersection point of a perpendicular of the line * through the point. */ Datum close_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); Point *result; - float8 invm; - LINE tmp; result = (Point *) palloc(sizeof(Point)); if (FPzero(line->B)) /* vertical? */ + point_construct(result, float8_div(-line->C, line->A), pt->y); + else if (FPzero(line->A)) /* horizontal? */ + point_construct(result, pt->x, float8_div(-line->C, line->B)); + else { - result->x = line->C; - result->y = pt->y; - PG_RETURN_POINT_P(result); - } - if (FPzero(line->A)) /* horizontal? */ - { - result->x = pt->x; - result->y = line->C; - PG_RETURN_POINT_P(result); - } - /* drop a perpendicular and find the intersection point */ + LINE perp; + + /* + * Drop a perpendicular and find the intersection point + * + * We need to invert and flip the sign on the slope to get the + * perpendicular. We might lose some precision on the division. It + * shouldn't be as much to turn the line. If it happens anyway, we + * will assume the point it on the line. + */ + line_construct_pm(&perp, pt, float8_div(line->B, line->A)); + if (!line_interpt_internal(result, &perp, line)) + *result = *pt; + } - /* invert and flip the sign on the slope to get a perpendicular */ - invm = float8_div(line->B, line->A); - line_construct_pm(&tmp, pt, invm); - line_interpt_internal(result, &tmp, line); PG_RETURN_POINT_P(result); } /* close_ps() * Closest point on line segment to specified point. * Take the closest endpoint if the point is left, right, * above, or below the segment, otherwise find the intersection * point of the segment and its perpendicular through the point. * -- 2.13.5 (Apple Git-94)