[postgis-tickets] r15896 - Fix division by zero in lw_dist2d_arc_arc

Paul Ramsey pramsey at cleverelephant.ca
Thu Oct 5 05:41:57 PDT 2017


Author: pramsey
Date: 2017-10-05 05:41:57 -0700 (Thu, 05 Oct 2017)
New Revision: 15896

Modified:
   trunk/liblwgeom/cunit/cu_measures.c
   trunk/liblwgeom/measures.c
Log:
Fix division by zero in lw_dist2d_arc_arc
(References #3879)



Modified: trunk/liblwgeom/cunit/cu_measures.c
===================================================================
--- trunk/liblwgeom/cunit/cu_measures.c	2017-10-05 11:33:02 UTC (rev 15895)
+++ trunk/liblwgeom/cunit/cu_measures.c	2017-10-05 12:41:57 UTC (rev 15896)
@@ -676,25 +676,130 @@
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
 
-	/* inscribed and closest on arcs */
+	/* Concentric: and fully parallel */
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
-	A1.x = -0.5; A1.y = 0.0;
-	A2.x =  0.0; A2.y = 0.5;
-	A3.x =  0.5; A3.y = 0.0;
+	A1.x = -2.0; A1.y = 0.0;
+	A2.x =  0.0; A2.y = 2.0;
+	A3.x =  2.0; A3.y = 0.0;
 	rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
-	//printf("distance %g\n", dl.distance);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
+	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
+
+	/* Concentric: with A fully included in B's range */
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -0.5 / sqrt(2.0); A1.y =  0.5 / sqrt(2.0);
+	A2.x =  0.0;             A2.y =  0.5;
+	A3.x =  0.5 / sqrt(2.0); A3.y =  0.5 / sqrt(2.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.5, 0.000001);
 
-	/* inscribed and closest not on arcs */
+	/* Concentric: with A partially included in B's range */
 	lw_dist2d_distpts_init(&dl, DIST_MIN);
-	A1.x = -0.5; A1.y =  0.0;
-	A2.x =  0.0; A2.y = -0.5;
-	A3.x =  0.5; A3.y =  0.0;
+	A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
+	A2.x = -0.5;             A2.y =  0.0;
+	A3.x = -0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
 	rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
-	//printf("distance %g\n", dl.distance);
 	CU_ASSERT_EQUAL( rv, LW_SUCCESS );
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
+
+	/* Concentric: with A and B without parallel segments */
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
+	A2.x =  0.0;             A2.y = -0.5;
+	A3.x =  0.5 / sqrt(2.0); A3.y = -0.5 / sqrt(2.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.736813, 0.000001);
+
+	/* Concentric: Arcs are the same */
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -1.0; A1.y =  0.0;
+	A2.x =  0.0; A2.y =  1.0;
+	A3.x =  1.0; A3.y =  0.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.0, 0.000001);
+
+	/* Check different orientations */
+	B1.x = -10.0;  B1.y =  0.0;
+	B2.x =   0.0 ; B2.y = 10.0;
+	B3.x =  10.0 ; B3.y =  0.0;
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -22.0; A1.y =   0.0;
+	A2.x = -17.0; A2.y =  -5.0;
+	A3.x = -12.0; A3.y =   0.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, 2.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -19.0; A1.y =   0.0;
+	A2.x = -14.0; A2.y =  -5.0;
+	A3.x = - 9.0; A3.y =   0.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, 1.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -9.0;  A1.y =   0.0;
+	A2.x = -4.0;  A2.y =  -5.0;
+	A3.x =  1.0;  A3.y =   0.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, 1.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -1.0;  A1.y =   0.0;
+	A2.x =  4.0;  A2.y =  -5.0;
+	A3.x =  9.0;  A3.y =   0.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, 1.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x =  1.0;  A1.y =   0.0;
+	A2.x =  6.0;  A2.y =  -5.0;
+	A3.x = 11.0;  A3.y =   0.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, 1.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = 11.0;  A1.y =   0.0;
+	A2.x = 16.0;  A2.y =  -5.0;
+	A3.x = 21.0;  A3.y =   0.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, 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;
+	A3.x = - 5.0; A3.y =  -6.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, 1.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -5.0; A1.y =  0.0;
+	A2.x =  0.0; A2.y =  5.0;
+	A3.x =  5.0; A3.y =  0.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, 5.0, 0.000001);
+
+	lw_dist2d_distpts_init(&dl, DIST_MIN);
+	A1.x = -5.0; A1.y =  0.0;
+	A2.x =  0.0; A2.y = -5.0;
+	A3.x =  5.0; A3.y =  0.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, 5.0, 0.000001);
+
+
 }
 
 static void

Modified: trunk/liblwgeom/measures.c
===================================================================
--- trunk/liblwgeom/measures.c	2017-10-05 11:33:02 UTC (rev 15895)
+++ trunk/liblwgeom/measures.c	2017-10-05 12:41:57 UTC (rev 15896)
@@ -1473,6 +1473,12 @@
 	return LW_TRUE;
 }
 
+/* Auxiliary function to calculate the distance between 2 concentric arcs*/
+int lw_dist2d_arc_arc_concentric(	const POINT2D *A1, const POINT2D *A2,
+					const POINT2D *A3, double radius_A,
+					const POINT2D *B1, const POINT2D *B2,
+					const POINT2D *B3, double radius_B,
+					const POINT2D *CENTER, DISTPTS *dl);
 
 int
 lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
@@ -1514,6 +1520,15 @@
 	if ( radius_B < 0 )
 		return lw_dist2d_seg_arc(B1, B3, A1, A2, A3, dl);
 
+	/* Center-center distance */
+	d = distance2d_pt_pt(&CA, &CB);
+
+	/* Concentric arcs */
+	if ( FP_EQUALS(d, 0.0) )
+		return lw_dist2d_arc_arc_concentric(A1, A2, A3, radius_A,
+						    B1, B2, B3, radius_B,
+						    &CA, dl);
+
 	/* Make sure that arc "A" has the bigger radius */
 	if ( radius_B > radius_A )
 	{
@@ -1525,15 +1540,6 @@
 		d = radius_B; radius_B = radius_A; radius_A = d;
 	}
 
-	/* Center-center distance */
-	d = distance2d_pt_pt(&CA, &CB);
-
-	/* Equal circles. Arcs may intersect at multiple points, or at none! */
-	if ( FP_EQUALS(d, 0.0) && FP_EQUALS(radius_A, radius_B) )
-	{
-		lwerror("lw_dist2d_arc_arc can't handle cojoint circles, uh oh");
-	}
-
 	/* Circles touch at a point. Is that point within the arcs? */
 	if ( d == (radius_A + radius_B) )
 	{
@@ -1655,6 +1661,143 @@
 	return LW_TRUE;
 }
 
+int
+lw_dist2d_arc_arc_concentric(	const POINT2D *A1, const POINT2D *A2,
+				const POINT2D *A3, double radius_A,
+				const POINT2D *B1, const POINT2D *B2,
+				const POINT2D *B3, double radius_B,
+				const POINT2D *CENTER, DISTPTS *dl)
+{
+	int seg_size;
+	double dist_sqr, shortest_sqr;
+	const POINT2D *P1;
+	const POINT2D *P2;
+	POINT2D proj;
+
+	if (radius_A == radius_B)
+	{
+		/* Check if B1 or B3 are in the same side as A2 in the A1-A3 arc */
+		seg_size = lw_segment_side(A1, A3, A2);
+		if (seg_size == lw_segment_side(A1, A3, B1))
+		{
+			dl->p1 = *B1;
+			dl->p2 = *B1;
+			dl->distance = 0;
+			return LW_TRUE;
+		}
+		if (seg_size == lw_segment_side(A1, A3, B3))
+		{
+			dl->p1 = *B3;
+			dl->p2 = *B3;
+			dl->distance = 0;
+			return LW_TRUE;
+		}
+		/* Check if A1 or A3 are in the same side as B2 in the B1-B3 arc */
+		seg_size = lw_segment_side(B1, B3, B2);
+		if (seg_size == lw_segment_side(B1, B3, A1))
+		{
+			dl->p1 = *A1;
+			dl->p2 = *A1;
+			dl->distance = 0;
+			return LW_TRUE;
+		}
+		if (seg_size == lw_segment_side(B1, B3, A3))
+		{
+			dl->p1 = *A3;
+			dl->p2 = *A3;
+			dl->distance = 0;
+			return LW_TRUE;
+		}
+	}
+	else
+	{
+		/* Check if any projection of B ends are in A*/
+		seg_size = lw_segment_side(A1, A3, A2);
+
+		/* B1 */
+		proj.x = CENTER->x + (B1->x - CENTER->x) * radius_A / radius_B;
+		proj.y = CENTER->y + (B1->y - CENTER->y) * radius_A / radius_B;
+
+		if (seg_size == lw_segment_side(A1, A3, &proj))
+		{
+			dl->p1 = proj;
+			dl->p2 = *B1;
+			dl->distance = fabs(radius_A - radius_B);
+			return LW_TRUE;
+		}
+		/* B3 */
+		proj.x = CENTER->x + (B3->x - CENTER->x) * radius_A / radius_B;
+		proj.y = CENTER->y + (B3->y - CENTER->y) * radius_A / radius_B;
+		if (seg_size == lw_segment_side(A1, A3, &proj))
+		{
+			dl->p1 = proj;
+			dl->p2 = *B3;
+			dl->distance = fabs(radius_A - radius_B);
+			return LW_TRUE;
+		}
+
+		/* Now check projections of A in B */
+		seg_size = lw_segment_side(B1, B3, B2);
+
+		/* A1 */
+		proj.x = CENTER->x + (A1->x - CENTER->x) * radius_B / radius_A;
+		proj.y = CENTER->y + (A1->y - CENTER->y) * radius_B / radius_A;
+		if (seg_size == lw_segment_side(B1, B3, &proj))
+		{
+			dl->p1 = proj;
+			dl->p2 = *A1;
+			dl->distance = fabs(radius_A - radius_B);
+			return LW_TRUE;
+		}
+
+		/* A3 */
+		proj.x = CENTER->x + (A3->x - CENTER->x) * radius_B / radius_A;
+		proj.y = CENTER->y + (A3->y - CENTER->y) * radius_B / radius_A;
+		if (seg_size == lw_segment_side(B1, B3, &proj))
+		{
+			dl->p1 = proj;
+			dl->p2 = *A3;
+			dl->distance = fabs(radius_A - radius_B);
+			return LW_TRUE;
+		}
+	}
+
+	/* Check the shortest between the distances of the 4 ends */
+	shortest_sqr = dist_sqr = distance2d_sqr_pt_pt(A1, B1);
+	P1 = A1;
+	P2 = B1;
+
+	dist_sqr = distance2d_sqr_pt_pt(A1, B3);
+	if (dist_sqr < shortest_sqr)
+	{
+		shortest_sqr = dist_sqr;
+		P1 = A1;
+		P2 = B3;
+	}
+
+	dist_sqr = distance2d_sqr_pt_pt(A3, B1);
+	if (dist_sqr < shortest_sqr)
+	{
+		shortest_sqr = dist_sqr;
+		P1 = A3;
+		P2 = B1;
+	}
+
+	dist_sqr = distance2d_sqr_pt_pt(A3, B3);
+	if (dist_sqr < shortest_sqr)
+	{
+		shortest_sqr = dist_sqr;
+		P1 = A3;
+		P2 = B3;
+	}
+
+	dl->p1 = *P1;
+	dl->p2 = *P2;
+	dl->distance = sqrt(shortest_sqr);
+
+	return LW_TRUE;
+}
+
 /**
 Finds the shortest distance between two segments.
 This function is changed so it is not doing any comparasion of distance



More information about the postgis-tickets mailing list