[postgis-tickets] r17182 - More robust geography distance

[Paul Ramsey]

Author: pramsey
Date: 2019-01-18 10:02:30 -0800 (Fri, 18 Jan 2019)
New Revision: 17182

Modified:
   branches/2.5/liblwgeom/lwgeodetic.c
Log:
More robust geography distance
References #4290


Modified: branches/2.5/liblwgeom/lwgeodetic.c
===================================================================
--- branches/2.5/liblwgeom/lwgeodetic.c	2019-01-18 17:25:55 UTC (rev 17181)
+++ branches/2.5/liblwgeom/lwgeodetic.c	2019-01-18 18:02:30 UTC (rev 17182)
@@ -36,6 +36,15 @@
 int gbox_geocentric_slow = LW_FALSE;
 
 /**
+* Utility function for ptarray_contains_point_sphere()
+*/
+static int
+point3d_equals(const POINT3D *p1, const POINT3D *p2)
+{
+	return FP_EQUALS(p1->x, p2->x) && FP_EQUALS(p1->y, p2->y) && FP_EQUALS(p1->z, p2->z);
+}
+
+/**
 * Convert a longitude to the range of -PI,PI
 */
 double longitude_radians_normalize(double lon)
@@ -3373,6 +3382,10 @@
 	POINT3D AC; /* Center point of A1/A2 */
 	double min_similarity, similarity;
 
+	/* Boundary case */
+	if (point3d_equals(A1, P) || point3d_equals(A2, P))
+		return LW_TRUE;
+
 	/* The normalized sum bisects the angle between start and end. */
 	vector_sum(A1, A2, &AC);
 	normalize(&AC);
@@ -3380,13 +3393,46 @@
 	/* The projection of start onto the center defines the minimum similarity */
 	min_similarity = dot_product(A1, &AC);
 
-	/* The projection of candidate p onto the center */
-	similarity = dot_product(P, &AC);
+	/* If the edge is sufficiently curved, use the dot product test */
+	if (fabs(1.0 - min_similarity) > 1e-10)
+	{
+		/* The projection of candidate p onto the center */
+		similarity = dot_product(P, &AC);
 
-	/* If the point is more similar than the end, the point is in the cone */
-	if ( similarity > min_similarity || fabs(similarity - min_similarity) < 2e-16 )
+		/* If the projection of the candidate is larger than */
+		/* the projection of the start point, the candidate */
+		/* must be closer to the center than the start, so */
+		/* therefor inside the cone */
+		if (similarity > min_similarity)
+		{
+			return LW_TRUE;
+		}
+		else
+		{
+			return LW_FALSE;
+		}
+	}
+	else
 	{
-		return LW_TRUE;
+		/* Where the edge is very narrow, the dot product test */
+		/* fails, but we can use the almost-planar nature of the */
+		/* problem space then to test if the vector from the */
+		/* candidate to the start point in a different direction */
+		/* to the vector from candidate to end point */
+		/* If so, then candidate is between start and end */
+		POINT3D PA1, PA2;
+		vector_difference(P, A1, &PA1);
+		vector_difference(P, A2, &PA2);
+		normalize(&PA1);
+		normalize(&PA2);
+		if (dot_product(&PA1, &PA2) < 0.0)
+		{
+			return LW_TRUE;
+		}
+		else
+		{
+			return LW_FALSE;
+		}
 	}
 	return LW_FALSE;
 }
@@ -3393,15 +3439,6 @@
 
 
 /**
-* Utility function for ptarray_contains_point_sphere()
-*/
-static int
-point3d_equals(const POINT3D *p1, const POINT3D *p2)
-{
-	return FP_EQUALS(p1->x, p2->x) && FP_EQUALS(p1->y, p2->y) && FP_EQUALS(p1->z, p2->z);
-}
-
-/**
 * Utility function for edge_intersects(), signum with a tolerance
 * in determining if the value is zero.
 */