[postgis-tickets] r17388 - Circular string distances have some failure modes, and test

Paul Ramsey pramsey at cleverelephant.ca
Fri Apr 12 03:21:18 PDT 2019


Author: pramsey
Date: 2019-04-12 15:21:18 -0700 (Fri, 12 Apr 2019)
New Revision: 17388

Modified:
   branches/2.3/NEWS
   branches/2.3/liblwgeom/cunit/cu_measures.c
   branches/2.3/liblwgeom/measures.c
Log:
Circular string distances have some failure modes, and test
suite should test everything bi-directionally.
References #4326


Modified: branches/2.3/NEWS
===================================================================
--- branches/2.3/NEWS	2019-04-12 22:16:43 UTC (rev 17387)
+++ branches/2.3/NEWS	2019-04-12 22:21:18 UTC (rev 17388)
@@ -4,7 +4,9 @@
   * Bug Fixes and Enchantments *
 
   - #4361, Fix postgis_type_name with (GEOMETRYM,3) (Matt Bretl)
+  - #4326, Fix circular arc distance calculation (Paul Ramsey)
 
+
 PostGIS 2.3.9
 2019/03/11
 

Modified: branches/2.3/liblwgeom/cunit/cu_measures.c
===================================================================
--- branches/2.3/liblwgeom/cunit/cu_measures.c	2019-04-12 22:16:43 UTC (rev 17387)
+++ branches/2.3/liblwgeom/cunit/cu_measures.c	2019-04-12 22:21:18 UTC (rev 17388)
@@ -31,7 +31,8 @@
 }
 
 #define DIST2DTEST(str1, str2, res) \
-	do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance2d_tolerance)
+	do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance2d_tolerance);\
+	do_test_mindistance_tolerance(str2, str1, res, __LINE__, lwgeom_mindistance2d_tolerance)
 #define DIST3DTEST(str1, str2, res) \
 	do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance3d_tolerance)
 
@@ -504,7 +505,7 @@
 	LWGEOM *out = NULL;
 	double measure = 105.0;
 	char *str;
-	
+
 	/* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
 	geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
 	out = lwgeom_locate_along(geom, measure, 0.0);
@@ -513,7 +514,7 @@
 	lwgeom_free(out);
 	CU_ASSERT_STRING_EQUAL("MULTIPOINT M (55.226131 55.226131 105)", str);
 	lwfree(str);
-	
+
 	/* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
 	geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 3 4 2, 9 4 3), (1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
 	out = lwgeom_locate_along(geom, measure, 0.0);
@@ -532,12 +533,12 @@
 	POINT2D P, A1, A2, A3;
 	int rv;
 
-	
+
 	/* Point within unit semicircle, 0.5 units from arc */
 	A1.x = -1; A1.y = 0;
 	A2.x = 0 ; A2.y = 1;
 	A3.x = 1 ; A3.y = 0;
-	P.x  = 0 ; P.y  = 0.5;	
+	P.x  = 0 ; P.y  = 0.5;
 
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
@@ -545,7 +546,7 @@
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
 
 	/* Point outside unit semicircle, 0.5 units from arc */
-	P.x  = 0 ; P.y  = 1.5;	
+	P.x  = 0 ; P.y  = 1.5;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
@@ -552,7 +553,7 @@
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
 
 	/* Point outside unit semicircle, sqrt(2) units from arc end point*/
-	P.x  = 0 ; P.y  = -1;	
+	P.x  = 0 ; P.y  = -1;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
@@ -559,7 +560,7 @@
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
 
 	/* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
-	P.x  = 1 ; P.y  = 1;	
+	P.x  = 1 ; P.y  = 1;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
@@ -566,7 +567,7 @@
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
 
 	/* Point on unit semicircle midpoint */
-	P.x  = 0 ; P.y  = 1;	
+	P.x  = 0 ; P.y  = 1;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
@@ -573,7 +574,7 @@
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
 
 	/* Point on unit semicircle endpoint */
-	P.x  = 1 ; P.y  = 0;	
+	P.x  = 1 ; P.y  = 0;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
@@ -594,7 +595,7 @@
 	rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
 	//printf("distance %g\n", dl.distance);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
-	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);	
+	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
 }
 
 static void
@@ -605,7 +606,7 @@
 	DISTPTS dl;
 	POINT2D A1, A2, B1, B2, B3;
 	int rv;
-	
+
 	/* Unit semicircle */
 	B1.x = -1; B1.y = 0;
 	B2.x = 0 ; B2.y = 1;
@@ -686,7 +687,7 @@
 	DISTPTS dl;
 	POINT2D A1, A2, A3, B1, B2, B3;
 	int rv;
-	
+
 	/* Unit semicircle at 0,0 */
 	B1.x = -1; B1.y = 0;
 	B2.x = 0 ; B2.y = 1;
@@ -826,7 +827,6 @@
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
 
-
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
 	A1.x = -15.0; A1.y =  -6.0;
 	A2.x = -10.0; A2.y =  -1.0;
@@ -851,6 +851,17 @@
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
 
+	/* Ticket #4326 */
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -1.0; A1.y =  4.0;
+	A2.x =  0.0; A2.y =  5.0;
+	A3.x =  1.0; A3.y =  4.0;
+	B1.x =  1.0; B1.y =  6.0;
+	B2.x =  6.0; B2.y =  1.0;
+	B3.x =  9.0; B3.y =  7.0;
+	rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
+	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
+	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0475666, 0.000001);
 
 }
 
@@ -861,7 +872,7 @@
 
 	POINT2D A1, A2, A3;
 	double d;
-	
+
 	/* Unit semicircle at 0,0 */
 	A1.x = -1; A1.y = 0;
 	A2.x = 0 ; A2.y = 1;
@@ -907,7 +918,7 @@
 
 	/* Unit semi-circle above X axis */
 	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
-	
+
 	/* Point at origin */
 	P.x = P.y = 0;
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
@@ -990,7 +1001,7 @@
 
 	/* Unit semi-circle above X axis */
 	lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
-	
+
 	/* Line above top of semi-circle */
 	lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
@@ -1007,7 +1018,7 @@
 	CU_ASSERT_STRING_EQUAL("lw_dist2d_ptarray_ptarrayarc called with non-arc input", cu_error_msg);
 
 	lwline_free(lwline2);
-	
+
 	/* Line along side of semi-circle */
 	lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
 	lw_dist2d_distpts_init(&dl, DIST_MIN);

Modified: branches/2.3/liblwgeom/measures.c
===================================================================
--- branches/2.3/liblwgeom/measures.c	2019-04-12 22:16:43 UTC (rev 17387)
+++ branches/2.3/liblwgeom/measures.c	2019-04-12 22:21:18 UTC (rev 17388)
@@ -1496,7 +1496,6 @@
 {
 	POINT2D CA, CB; /* Center points of arcs A and B */
 	double radius_A, radius_B, d; /* Radii of arcs A and B */
-	POINT2D P; /* Temporary point P */
 	POINT2D D; /* Mid-point between the centers CA and CB */
 	int pt_in_arc_A, pt_in_arc_B; /* Test whether potential intersection point is within the arc */
 
@@ -1542,11 +1541,13 @@
 	if ( radius_B > radius_A )
 	{
 		const POINT2D *tmp;
+		POINT2D TP; /* Temporary point P */
+		double td;
 		tmp = B1; B1 = A1; A1 = tmp;
 		tmp = B2; B2 = A2; A2 = tmp;
 		tmp = B3; B3 = A3; A3 = tmp;
-		P = CB; CB = CA; CA = P;
-		d = radius_B; radius_B = radius_A; radius_A = d;
+		TP = CB; CB = CA; CA = TP;
+		td = radius_B; radius_B = radius_A; radius_A = td;
 	}
 
 	/* Circles touch at a point. Is that point within the arcs? */
@@ -1662,8 +1663,8 @@
 	{
 		lw_dist2d_pt_pt(A1, B1, dl);
 		lw_dist2d_pt_pt(A1, B3, dl);
-		lw_dist2d_pt_pt(A2, B1, dl);
-		lw_dist2d_pt_pt(A2, B3, dl);
+		lw_dist2d_pt_pt(A3, B1, dl);
+		lw_dist2d_pt_pt(A3, B3, dl);
 		return LW_TRUE;
 	}
 



More information about the postgis-tickets mailing list