From dd5479d85de0d5291c0b79606ef3b3571f5f4d44 Mon Sep 17 00:00:00 2001
From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Wed, 5 Feb 2014 23:00:07 +0000
Subject: [PATCH] #2534, st_distance is returning incorrect results for large
 geographies

git-svn-id: http://svn.osgeo.org/postgis/trunk@12228 b70326c6-7e19-0410-871a-916f4a2858ee
---
 NEWS                          |  1 +
 liblwgeom/cunit/cu_geodetic.c | 25 +++++++++++++++++++++++++
 liblwgeom/lwgeodetic.c        | 19 ++++++++++++++++---
 3 files changed, 42 insertions(+), 3 deletions(-)

diff --git a/NEWS b/NEWS
index e97a1244c..7ffde6541 100644
--- a/NEWS
+++ b/NEWS
@@ -43,6 +43,7 @@ PostGIS 2.2.0
   - #2512, Support for foreign tables and materialized views in
            raster_columns and raster_overviews
   - #2532, Add missing raster/geometry commutator operators
+  - #2534, st_distance is returning incorrect results for large geographies
   - #2540, Change GUC name for GDAL_DATA to postgis.gdal_datapath
   - #2543, invalid join selectivity error from simple query
   - #2552, Fix handling of NULL raster in ST_AsTIFF, ST_AsPNG
diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 2a22b1ff4..f9b41d542 100644
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -674,6 +674,31 @@ static void test_edge_intersects(void)
 	line2pts("LINESTRING(90.0 80.0, -90.0 90.0)", &B1, &B2);
 	rv = edge_intersects(&A1, &A2, &B1, &B2);
 	CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_LEFT|PIR_A_TOUCH_RIGHT) );
+
+	/* Antipodal straddles. Great circles cross but at opposite */
+	/* sides of the globe */
+	/* #2534 */
+	/* http://www.gcmap.com/mapui?P=60N+90E-20S+90E%0D%0A0N+0E-90.04868865037885W+57.44011727050777S%0D%0A&MS=wls&DU=mi */
+	line2pts("LINESTRING(90.0 60.0, 90.0 -20.0)", &A1, &A2);
+	line2pts("LINESTRING(0.0 0.0, -90.04868865037885 -57.44011727050777)", &B1, &B2);
+	rv = edge_intersects(&A1, &A2, &B1, &B2);
+	CU_ASSERT(rv == 0);
+
+	line2pts("LINESTRING(-5 0, 5 0)", &A1, &A2);
+	line2pts("LINESTRING(179 -5, 179 5)", &B1, &B2);
+	rv = edge_intersects(&A1, &A2, &B1, &B2);
+	CU_ASSERT(rv == 0);
+
+	line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
+	line2pts("LINESTRING(65 0, -105 0)", &B1, &B2);
+	rv = edge_intersects(&A1, &A2, &B1, &B2);
+	CU_ASSERT(rv == 0);
+	
+	line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
+	line2pts("LINESTRING(45 0, -125 0)", &B1, &B2);
+	rv = edge_intersects(&A1, &A2, &B1, &B2);
+	CU_ASSERT(rv == 0);
+	
 }
 
 static void test_edge_distance_to_point(void)
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index a542a870e..bc8bf0594 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -3078,7 +3078,7 @@ dot_product_side(const POINT3D *p, const POINT3D *q)
 int 
 edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
 {
-	POINT3D AN, BN;  /* Normals to plane A and plane B */
+	POINT3D AN, BN, VN;  /* Normals to plane A and plane B */
 	double ab_dot;
 	int a1_side, a2_side, b1_side, b2_side;
 	int rv = PIR_NO_INTERACT;
@@ -3126,8 +3126,21 @@ edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const P
 	if ( a1_side != a2_side && (a1_side + a2_side) == 0 &&
 	     b1_side != b2_side && (b1_side + b2_side) == 0 )
 	{
-		/* Mid-point intersection! */
-		return PIR_INTERSECTS;
+		/* Have to check if intersection point is inside both arcs */
+		unit_normal(&AN, &BN, &VN);
+		if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
+		{
+			return PIR_INTERSECTS;
+		}
+
+		/* Have to check if intersection point is inside both arcs */
+		vector_scale(&VN, -1);
+		if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
+		{
+			return PIR_INTERSECTS;
+		}
+		
+		return PIR_NO_INTERACT;
 	}
 
 	/* The rest are all intersects variants... */
-- 
2.40.0