From: Paul Ramsey Date: Thu, 5 Oct 2017 13:00:27 +0000 (+0000) Subject: division by zero in lw_dist2d_arc_arc X-Git-Tag: 2.4.1~19 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=41084cbc94189096f9e9d5259441ca2482c4bf29;p=postgis division by zero in lw_dist2d_arc_arc References #3879 git-svn-id: http://svn.osgeo.org/postgis/branches/2.4@15897 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index 89579838a..98b5fcfaa 100644 --- a/NEWS +++ b/NEWS @@ -6,6 +6,7 @@ YYYY/MM/DD - #3869, Fix build with "gold" linker - #3845, Gracefully handle short-measure issue - #3871, Performance tweak for geometry cmp function + - #3879, Division by zero in some arc cases PostGIS 2.4.0 2017/09/30 diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c index 08857abd3..dc714585e 100644 --- a/liblwgeom/cunit/cu_measures.c +++ b/liblwgeom/cunit/cu_measures.c @@ -676,25 +676,130 @@ test_lw_dist2d_arc_arc(void) 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); + 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); - //printf("distance %g\n", dl.distance); 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 diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index d36e9e8dd..d6fe302f8 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -1473,6 +1473,12 @@ lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const P 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 @@ lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, 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 @@ lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, 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 @@ lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, 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