0004-line-fixes-v06.patch

From e93b101fc83114a37b5a5b0f9cbd866d51a6820b Mon Sep 17 00:00:00 2001
From: Emre Hasegeli <emre@hasegeli.com>
Date: Sun, 28 May 2017 11:35:17 +0200
Subject: [PATCH 4/6] line-fixes-v06

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
* 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 | 231 +++++++++++++++++++++++++++-------------
 1 file changed, 155 insertions(+), 76 deletions(-)

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 16f853ae9b..a33ff3f56b 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -47,20 +47,22 @@ static int	point_inside(Point *p, int npts, Point *plist);
 
 /* Routines for two-dimensional lines */
 static inline void line_construct_pm(LINE *result, Point *pt, float8 m);
 static inline void line_construct_pts(LINE *result, Point *pt1, Point *pt2);
 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 two-dimensional line segments */
 static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
+static bool lseg_parallel_line(LSEG *lseg, LINE *line);
+static bool lseg_parallel_lseg(LSEG *l1, LSEG *l2);
 static bool	lseg_interpt_line(Point *result, LSEG *lseg, LINE *line);
 static bool	lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2);
 static int	lseg_crossing(float8 x, float8 y, float8 px, float8 py);
 static bool	lseg_contain_point(LSEG *lseg, Point *point);
 static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt);
 static float8 lseg_closept_line(Point *result, LSEG *lseg, LINE *line);
 static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2);
 
 /* Routines for two-dimensional boxes */
 static inline void box_construct(BOX *result, Point *pt1, Point *pt2);
@@ -970,20 +972,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);
@@ -1064,45 +1071,59 @@ line_construct_pp(PG_FUNCTION_ARGS)
 
 Datum
 line_intersect(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
 
 	PG_RETURN_BOOL(line_interpt_line(NULL, l1, l2));
 }
 
+
+/*
+ * Check whether the two lines are parallel
+ *
+ * They are parallel if their slopes are FPeq().
+ */
 Datum
 line_parallel(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
+	float8		m1,
+				m2;
 
-	if (FPzero(l1->B))
-		PG_RETURN_BOOL(FPzero(l2->B));
+	m1 = slope(0.0, l1->B, l1->A, 0.0);
+	m2 = slope(0.0, l2->B, l2->A, 0.0);
 
-	PG_RETURN_BOOL(FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))));
+	PG_RETURN_BOOL(FPeq(m1, m2));
 }
 
+
+/*
+ * Check whether the two lines are perpendicular
+ *
+ * They are perpendicular if the slope of one is FPeq() with the inverse
+ * slope of the other.
+ */
 Datum
 line_perp(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
+	float8		m1,
+				m2;
 
-	if (FPzero(l1->A))
-		PG_RETURN_BOOL(FPzero(l2->B));
-	else if (FPzero(l1->B))
-		PG_RETURN_BOOL(FPzero(l2->A));
+	m1 = slope(0.0, l1->B, l1->A, 0.0);
+	m2 = slope(l2->A, 0.0, l2->B, 0.0);
 
-	PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B),
-								   float8_mul(l1->B, l2->A)), -1.0));
+	PG_RETURN_BOOL(FPeq(m1, m2));
 }
 
 Datum
 line_vertical(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->B));
 }
 
@@ -1114,25 +1135,28 @@ 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;
 
-	if (!FPzero(l2->A) && !isnan(l2->A))
+	if (!FPzero(l1->A) && !isnan(l1->A) &&
+		!FPzero(l2->A) && !isnan(l2->A))
 		ratio = float8_div(l1->A, l2->A);
-	else if (!FPzero(l2->B) && !isnan(l2->B))
+	else if (!FPzero(l1->B) && !isnan(l1->B) &&
+			 !FPzero(l2->B) && !isnan(l2->B))
 		ratio = float8_div(l1->B, l2->B);
-	else if (!FPzero(l2->C) && !isnan(l2->C))
+	else if (!FPzero(l1->C) && !isnan(l1->C) &&
+			 !FPzero(l2->C) && !isnan(l2->C))
 		ratio = float8_div(l1->C, l2->C);
 	else
 		ratio = 1.0;
 
 	PG_RETURN_BOOL(((isnan(l1->A) && isnan(l2->A)) ||
 					FPeq(l1->A, float8_mul(ratio, l2->A))) &&
 				   ((isnan(l1->B) && isnan(l2->B)) ||
 					FPeq(l1->B, float8_mul(ratio, l2->B))) &&
 				   ((isnan(l1->C) && isnan(l2->C)) ||
 					FPeq(l1->C, float8_mul(ratio, l2->C))));
@@ -1144,30 +1168,35 @@ 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_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);
@@ -1194,40 +1223,56 @@ line_interpt(PG_FUNCTION_ARGS)
  * If the lines have NaN constants, we will return true, and the intersection
  * point would have NaN coordinates.  We shouldn't return false in this case
  * because that would mean the lines are parallel.
  */
 static bool
 line_interpt_line(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(l2->A))		/* l2 horizontal? */
+	{
+		y = float8_div(-l2->C, l2->B);
+		x = float8_div(-float8_pl(float8_mul(l1->B, y), l1->C), l1->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->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
 	{
-		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;
 
-		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;
 }
 
 
 /***********************************************************************
  **
@@ -2021,59 +2066,83 @@ lseg_length(PG_FUNCTION_ARGS)
 Datum
 lseg_intersect(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 
 	PG_RETURN_BOOL(lseg_interpt_lseg(NULL, l1, l2));
 }
 
 
+/*
+ * Check whether the line segment is parallel with the line
+ *
+ * They are parallel if their slopes are FPeq().
+ */
+static bool
+lseg_parallel_line(LSEG *lseg, LINE *line)
+{
+	float8		m1,
+				m2;
+
+	m1 = slope(lseg->p[0].x, lseg->p[1].x, lseg->p[0].y, lseg->p[1].y);
+	m2 = slope(0.0, line->B, line->A, 0.0);
+
+	return FPeq(m1, m2);
+}
+
+
+/*
+ * Check whether the two line segment are parallel
+ *
+ * They are parallel if their slopes are FPeq().
+ */
+static bool
+lseg_parallel_lseg(LSEG *l1, LSEG *l2)
+{
+	float8		m1,
+				m2;
+
+	m1 = slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y);
+	m2 = slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y);
+
+	return FPeq(m1, m2);
+}
+
 Datum
 lseg_parallel(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 
-	PG_RETURN_BOOL(FPeq(slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y),
-						slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y)));
+	PG_RETURN_BOOL(lseg_parallel_lseg(l1, l2));
 }
 
+
 /* 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 slope() and the results seem better.
- * - thomas 1998-01-31
+ * They are perpendicular if the slope of one is FPeq() with the inverse
+ * slope of the other.
  */
 Datum
 lseg_perp(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 	float8		m1,
 				m2;
 
 	m1 = slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y);
-	m2 = slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y);
+	m2 = 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(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));
 }
 
@@ -2321,30 +2390,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);
@@ -2509,36 +2570,40 @@ 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;
 	float8		m;
+	Point		closept;
 	LINE		tmp;
 
 	/* Drop a perpendicular and find the intersection point */
 	m = slope(line->A, 0.0, line->B, 0.0);
 	line_construct_pm(&tmp, point, m);
-	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))
+		closept = *point;
+
+	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));
@@ -2590,64 +2655,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_parallel_lseg(l1, 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);
 }
 
 
 /*
@@ -2808,20 +2884,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_parallel_line(lseg, 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.13.6 (Apple Git-96)