From 83049a69b87fc6e800e00823eb4424bcaa9bee90 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-v04 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 to line operators crashing on floating point precision loss The changes are also aiming make line operators more symmetric and less sensitive to floating point precision loss. The EPSILON interferes with every minor change in different ways. It is hard to guess which behaviour is more expected, but we can assume threating both sides of the operators more equally is more 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 | 108 ++++++++++++++++++++++++++-------------- 1 file changed, 70 insertions(+), 38 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index a4f6c1992e..7a444ad7a5 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 (FPzero(line->A) && FPzero(line->B)) + 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)); } @@ -2667,34 +2695,38 @@ lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line) Datum close_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); Point *result; result = (Point *) palloc(sizeof(Point)); if (FPzero(line->B)) /* vertical? */ - point_construct(result, line->C, pt->y); + point_construct(result, float8_div(-line->C, line->A), pt->y); else if (FPzero(line->A)) /* horizontal? */ - point_construct(result, pt->x, line->C); + point_construct(result, pt->x, float8_div(-line->C, line->B)); else { LINE tmp; /* * Drop a perpendicular and find the intersection point * - * We need to invert the slope to get the perpendicular. + * We need to invert 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 + * is on the line. */ line_construct_pm(&tmp, pt, float8_div(line->B, line->A)); - line_interpt_internal(result, &tmp, line); + if (!line_interpt_internal(result, &tmp, line)) + *result = *pt; } 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 -- 2.13.5 (Apple Git-94)