[postgis-tickets] r15899 - division by zero in lw_dist2d_arc_arc
Paul Ramsey
pramsey at cleverelephant.ca
Thu Oct 5 06:04:28 PDT 2017
Author: pramsey
Date: 2017-10-05 06:04:27 -0700 (Thu, 05 Oct 2017)
New Revision: 15899
Modified:
branches/2.2/NEWS
branches/2.2/liblwgeom/cunit/cu_measures.c
branches/2.2/liblwgeom/measures.c
Log:
division by zero in lw_dist2d_arc_arc
Closes #3879
Modified: branches/2.2/NEWS
===================================================================
--- branches/2.2/NEWS 2017-10-05 13:00:35 UTC (rev 15898)
+++ branches/2.2/NEWS 2017-10-05 13:04:27 UTC (rev 15899)
@@ -9,6 +9,7 @@
not a polygon with ST_ConcaveHull
- #3866, Rare crash generating TWKB with large coordinate values
- #3869, Fix build with "gold" linker
+ - #3879, Division by zero in some arc cases
PostGIS 2.2.5
Modified: branches/2.2/liblwgeom/cunit/cu_measures.c
===================================================================
--- branches/2.2/liblwgeom/cunit/cu_measures.c 2017-10-05 13:00:35 UTC (rev 15898)
+++ branches/2.2/liblwgeom/cunit/cu_measures.c 2017-10-05 13:04:27 UTC (rev 15899)
@@ -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: branches/2.2/liblwgeom/measures.c
===================================================================
--- branches/2.2/liblwgeom/measures.c 2017-10-05 13:00:35 UTC (rev 15898)
+++ branches/2.2/liblwgeom/measures.c 2017-10-05 13:04:27 UTC (rev 15899)
@@ -1459,6 +1459,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,
@@ -1500,6 +1506,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 )
{
@@ -1511,15 +1526,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) )
{
@@ -1641,6 +1647,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