From 7c0daa810e2d4500638d2c7062c51f1dc9dcefcd 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-v05

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
* Don't return NULL for closest point to line segment
* Fix closest point to line segment from line 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 | 115 ++++++++++++++++++++++++----------------
 1 file changed, 68 insertions(+), 47 deletions(-)

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 2daa8e12cb..77dc0df5a7 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -972,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);
@@ -1090,25 +1095,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));
 }
 
@@ -1120,25 +1128,25 @@ 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))
+	if (!FPzero(l1->A) && !FPzero(l2->A))
 		ratio = float8_div(l1->A, l2->A);
-	else if (!FPzero(l2->B))
+	else if (!FPzero(l1->B) && !FPzero(l2->B))
 		ratio = float8_div(l1->B, l2->B);
-	else if (!FPzero(l2->C))
+	else if (!FPzero(l1->C) && !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)) &&
 				   FPeq(l1->B, float8_mul(ratio, l2->B)) &&
 				   FPeq(l1->C, float8_mul(ratio, l2->C)));
 }
 
 
@@ -1147,30 +1155,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_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);
+
+	if (!FPzero(l1->A) && !FPzero(l2->A))
+		ratio = float8_div(l1->A, l2->A);
+	else if (!FPzero(l1->B) && !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))),
+								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);
@@ -1193,38 +1206,56 @@ 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(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 = 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
 
@@ -2027,46 +2058,35 @@ lseg_parallel(PG_FUNCTION_ARGS)
 
 	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);
 
 	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 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 = 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(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));
 }
 
@@ -2627,34 +2647,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
@@ -2756,24 +2780,24 @@ close_ps(PG_FUNCTION_ARGS)
 	 */
 	line_construct_pm(&tmp, pt, invm);
 #ifdef GEODEBUG
 	printf("close_ps- tmp A %f  B %f   C %f\n",
 		   tmp->A, tmp->B, tmp->C);
 #endif
 
 	/*
 	 * 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 a SQL NULL if so.
+	 * roundoff issues.  In this case, assume the point is on the line segment.
 	 */
 	if (!lseg_interpt_line_internal(result, lseg, &tmp))
-		PG_RETURN_NULL();
+		*result = *pt;
 
 #ifdef GEODEBUG
 	printf("close_ps- result.x %f  result.y %f\n", result->x, result->y);
 #endif
 	PG_RETURN_POINT_P(result);
 }
 
 
 /* close_lseg()
  * Closest point to l1 on l2.
@@ -2799,28 +2823,25 @@ close_lseg(PG_FUNCTION_ARGS)
 	dists2[0] = dist_ps_internal(&l2->p[0], l1);
 	dists2[1] = dist_ps_internal(&l2->p[1], l1);
 	if (float8_le(dists2[0], dists2[1]))
 		closer_end2 = 0;
 	else
 		closer_end2 = 1;
 
 	if (float8_le(dists1[closer_end1], dists2[closer_end2]))
 		point = PointPGetDatum(&l1->p[closer_end1]);
 	else
-	{
 		point = DirectFunctionCall2(close_ps,
 									PointPGetDatum(&l2->p[closer_end2]),
 									LsegPGetDatum(l1));
-		point = DirectFunctionCall2(close_ps, point, LsegPGetDatum(l2));
-	}
 
-	return point;
+	return DirectFunctionCall2(close_ps, point, LsegPGetDatum(l2));
 }
 
 /* close_pb()
  * Closest point on or in box to specified point.
  */
 Datum
 close_pb(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
-- 
2.13.6 (Apple Git-96)

