From 00fad9622bd5d6ce1b0e20c8662149b6e2a81e7e Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sun, 28 May 2017 11:35:17 +0200 Subject: [PATCH 4/6] line-fixes-v07 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 * Reject same points on line_construct_pp() * 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 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 | 143 +++++++++++++++++++++++++++------------- 1 file changed, 98 insertions(+), 45 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index c925a99876..24d2602f4b 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -40,20 +40,21 @@ static inline void point_add_point(Point *result, Point *pt1, Point *pt2); static inline void point_sub_point(Point *result, Point *pt1, Point *pt2); static inline void point_mul_point(Point *result, Point *pt1, Point *pt2); static inline void point_div_point(Point *result, Point *pt1, Point *pt2); static inline bool point_eq_point(Point *pt1, Point *pt2); static inline float8 point_dt(Point *pt1, Point *pt2); static inline float8 point_sl(Point *pt1, Point *pt2); 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); static float8 line_closept_point(Point *result, LINE *line, Point *pt); /* Routines for line segments */ static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2); static inline float8 lseg_sl(LSEG *lseg); static inline float8 lseg_invsl(LSEG *lseg); static bool lseg_interpt_line(Point *result, LSEG *lseg, LINE *line); @@ -972,20 +973,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_BINARY_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); @@ -1036,20 +1042,25 @@ line_construct(LINE *result, Point *pt, float8 m) /* line_construct_pp() * two points */ Datum line_construct_pp(PG_FUNCTION_ARGS) { Point *pt1 = PG_GETARG_POINT_P(0); 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); } /*---------------------------------------------------------- * Relative position routines. *---------------------------------------------------------*/ @@ -1072,22 +1083,26 @@ 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(l2->A)) + PG_RETURN_BOOL(FPzero(l1->B)); else if (FPzero(l1->B)) PG_RETURN_BOOL(FPzero(l2->A)); + else 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)); } Datum line_vertical(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); @@ -1131,20 +1146,34 @@ line_eq(PG_FUNCTION_ARGS) FPeq(l1->B, float8_mul(ratio, l2->B))) && ((isnan(l1->C) && isnan(l2->C)) || FPeq(l1->C, float8_mul(ratio, l2->C)))); } /*---------------------------------------------------------- * 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 */ static inline float8 line_invsl(LINE *line) { if (FPzero(line->A)) return DBL_MAX; if (FPzero(line->B)) return 0.0; @@ -1153,30 +1182,35 @@ line_invsl(LINE *line) /* 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_line(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 = 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() * 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); @@ -1208,35 +1242,39 @@ static bool line_interpt_line(Point *result, LINE *l1, LINE *l2) { float8 x, y; 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->B)) /* l2 vertical? */ { - x = l2->C; - y = float8_pl(float8_mul(l1->A, x), l1->C); + x = float8_div(-l2->C, l2->A); + y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B); } else { + /* XXX This check break commutative property. */ if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->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(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); return true; } /*********************************************************************** ** @@ -2345,30 +2383,22 @@ dist_pb(PG_FUNCTION_ARGS) } /* * Distance from a lseg to a line */ 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(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)); } /* * Distance from a lseg to a box */ Datum dist_sb(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); BOX *box = PG_GETARG_BOX_P(1); @@ -2533,34 +2563,43 @@ lseg_interpt_line(Point *result, LSEG *lseg, LINE *line) /* * The intersection point of a perpendicular of the line * through the point. * * This sets the closest point to the *result if it is not NULL and returns * the distance to the closest point. */ static float8 line_closept_point(Point *result, LINE *line, Point *point) { - bool retval; + Point closept; LINE tmp; /* Drop a perpendicular and find the intersection point */ line_construct(&tmp, point, line_invsl(line)); - retval = line_interpt_line(result, line, &tmp); - Assert(retval); /* XXX We need something better. */ /* - * XXX We could use the distance to the closest point, but - * line_interpt_line() is currently giving wrong results. + * Ordinarily we should always find an intersection point, but that could + * fail in the presence of NaN coordinates, and perhaps even from simple + * roundoff issues. */ - return fabs((line->A * point->x + line->B * point->y + line->C) / - HYPOT(line->A, line->B)); + if (!line_interpt_line(&closept, &tmp, line)) + { + if (result != NULL) + *result = *point; + + return get_float8_nan(); + } + + if (result != NULL) + *result = closept; + + return point_dt(&closept, point); } 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)); @@ -2610,64 +2649,75 @@ close_ps(PG_FUNCTION_ARGS) PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } /* * Closest point to l1 on l2. * * This sets the closest point to the *result if it is not NULL and returns - * the distance to the closest point. - * - * XXX This function is wrong. If must never set the *result to a point on - * the second segment. + * the distance to the closest point. We first eliminate the case + * the segments intersecting with each other. Then we try to find + * the closest point on the first segment by trying the endpoints of + * the second. Though, it is still possible for the closest point to be + * one of the endpoints, so we test them. */ static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2) { Point point; float8 dist, d; - d = lseg_closept_point(NULL, l1, &l2->p[0]); - dist = d; - if (result != NULL) - *result = l2->p[0]; + if (lseg_interpt_lseg(result, l1, l2)) + return 0.0; - d = lseg_closept_point(NULL, l1, &l2->p[1]); + 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); - + 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; } Datum close_lseg(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); 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))) PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } /* @@ -2828,20 +2878,23 @@ lseg_closept_line(Point *result, LSEG *lseg, LINE *line) } } Datum close_ls(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); 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))) PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } /* -- 2.14.3 (Apple Git-98)