From 0f8e84f7576a6ff0e4b3c9abd6edba9a79490039 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-v10

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()
* Fix perpendicular operator when negative values are involved
* 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.

Discussion: https://www.postgresql.org/message-id/CAE2gYzxF7-5djV6-cEvqQu-fNsnt%3DEqbOURx7ZDg%2BVv6ZMTWbg%40mail.gmail.com

Parallel discussions:

* https://www.postgresql.org/message-id/CAE2gYzw_-z%3DV2kh8QqFjenu%3D8MJXzOP44wRW%3DAzzeamrmTT1%3DQ%40mail.gmail.com
* https://www.postgresql.org/message-id/20180201.205138.34583581.horiguchi.kyotaro@lab.ntt.co.jp
---
 src/backend/utils/adt/geo_ops.c | 159 ++++++++++++++++++++++++++--------------
 1 file changed, 106 insertions(+), 53 deletions(-)

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index df185affcb..d5d4797b3f 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);
@@ -1029,20 +1035,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.
  *---------------------------------------------------------*/
 
@@ -1065,25 +1076,29 @@ 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(l2->A))
+		PG_RETURN_BOOL(FPzero(l1->B));
+	if (FPzero(l1->B))
 		PG_RETURN_BOOL(FPzero(l2->A));
+	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_div(float8_mul(l1->A, l2->A),
+								   float8_mul(l1->B, l2->B)), -1.0));
 }
 
 Datum
 line_vertical(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->B));
 }
 
@@ -1124,20 +1139,34 @@ line_eq(PG_FUNCTION_ARGS)
 				   (float8_eq(l1->A, l2->A) &&
 					float8_eq(l1->B, l2->B) &&
 					float8_eq(l1->C, 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;
@@ -1146,30 +1175,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))	/* intersecting? */
 		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);
@@ -1196,41 +1230,44 @@ 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(l2->B))		/* l2 vertical? */
-			return false;
-
-		x = l1->C;
-		y = float8_pl(float8_mul(l2->A, x), l2->C);
-	}
-	else if (FPzero(l2->B))		/* l2 vertical? */
-	{
-		x = l2->C;
-		y = float8_pl(float8_mul(l1->A, x), l1->C);
-	}
-	else
+	if (!FPzero(l1->B))
 	{
 		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);
 	}
+	else if (!FPzero(l2->B))
+	{
+		if (FPeq(l1->A, float8_mul(l2->A, float8_div(l1->B, l2->B))))
+			return false;
+
+		x = float8_div(float8_mi(float8_mul(l2->B, l1->C),
+								 float8_mul(l1->B, l2->C)),
+					   float8_mi(float8_mul(l2->A, l1->B),
+								 float8_mul(l1->A, l2->B)));
+		y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B);
+	}
+	else
+		return false;
 
 	if (result != NULL)
 		point_construct(result, x, y);
 
 	return true;
 }
 
 
 /***********************************************************************
  **
@@ -2336,30 +2373,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,33 +2562,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;
-	Point       closept;
+	Point		closept;
 	LINE		tmp;
 
-	/* We drop a perpendicular to find the intersection point. */
+	/*
+	 * We drop a perpendicular to find the intersection point.  Ordinarily
+	 * we should always find it, but that can fail in the presence of NaN
+	 * coordinates, and perhaps even from simple roundoff issues.
+	 */
 	line_construct(&tmp, point, line_invsl(line));
-	retval = line_interpt_line(&closept, line, &tmp);
-	Assert(retval);		/* XXX: We need something better. */
+	if (!line_interpt_line(&closept, &tmp, line))
+	{
+		if (result != NULL)
+			*result = *point;
+
+		return get_float8_nan();
+	}
 
 	if (result != NULL)
 		*result = closept;
 
-	/* Then we calculate the distance between the points. */
 	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;
 
@@ -2619,52 +2655,66 @@ close_ps(PG_FUNCTION_ARGS)
  * This sets the closest point to the *result if it is not NULL and returns
  * the distance to the closest point.
  */
 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];
+	/* First, we handle the case when the line segments are intersecting. */
+	if (lseg_interpt_lseg(result, l1, l2))
+		return 0.0;
 
-	d = lseg_closept_point(NULL, l1, &l2->p[1]);
+	/*
+	 * Then, we find the closest points from the endpoints of the second
+	 * line segment, and keep the closest one.
+	 */
+	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);
-
+	/* The closest point can still be one of the endpoints, so we test them. */
+	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);
 }
 
 
 /*
@@ -2825,20 +2875,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)

