From: Paul Ramsey Date: Thu, 27 Sep 2012 20:23:03 +0000 (+0000) Subject: Distance calculation support for arc features (#2018). Commit adds lowest level primi... X-Git-Tag: 2.1.0beta2~616 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=503d866d8ac13b3b9d3e5a6b300c5f770caddbe6;p=postgis Distance calculation support for arc features (#2018). Commit adds lowest level primitive support for distance calculations on single arcs. git-svn-id: http://svn.osgeo.org/postgis/trunk@10335 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c index 7f53c20d6..963dcfb2a 100644 --- a/liblwgeom/cunit/cu_measures.c +++ b/liblwgeom/cunit/cu_measures.c @@ -345,6 +345,203 @@ test_lwgeom_locate_along(void) lwfree(str); } +static void +test_lw_dist2d_pt_arc(void) +{ + /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */ + DISTPTS dl; + 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; + + lw_dist2d_distpts_init(&dl, DIST_MIN); + rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl); + 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; + lw_dist2d_distpts_init(&dl, DIST_MIN); + rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl); + 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; + lw_dist2d_distpts_init(&dl, DIST_MIN); + rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl); + 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; + lw_dist2d_distpts_init(&dl, DIST_MIN); + rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001); + + /* Point on unit semicircle midpoint */ + 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_DOUBLE_EQUAL(dl.distance, 0, 0.000001); + + /* Point on unit semicircle endpoint */ + 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_DOUBLE_EQUAL(dl.distance, 0, 0.000001); + + /* Point inside closed circle */ + P.x = 0 ; P.y = 0.5; + A2.x = 1; A2.y = 0; + A3 = A1; + lw_dist2d_distpts_init(&dl, DIST_MIN); + rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl); + //printf("distance %g\n", dl.distance); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001); +} + +static void +test_lw_dist2d_seg_arc(void) +{ + /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */ + + DISTPTS dl; + POINT2D A1, A2, B1, B2, B3; + int rv; + + /* Unit semicircle */ + B1.x = -1; B1.y = 0; + B2.x = 0 ; B2.y = 1; + B3.x = 1 ; B3.y = 0; + + /* Edge above the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -2; A1.y = 2; + A2.x = 2 ; A2.y = 2; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001); + + /* Edge to the right of the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = 2; A1.y = -2; + A2.x = 2; A2.y = 2; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001); + + /* Edge to the left of the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -2; A1.y = -2; + A2.x = -2; A2.y = 2; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001); + + /* Edge within the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = 0; A1.y = 0; + A2.x = 0; A2.y = 0.5; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001); + + /* Edge grazing the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -2; A1.y = 1; + A2.x = 2; A2.y = 1; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001); + + /* Line grazing the unit semicircle, but edge not */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = 1; A1.y = 1; + A2.x = 2; A2.y = 1; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001); + + /* Edge intersecting the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = 0; A1.y = 0; + A2.x = 2; A2.y = 2; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001); + + /* Line intersecting the unit semicircle, but edge not */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -1; A1.y = 1; + A2.x = -2; A2.y = 2; + rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl); + //printf("distance %g\n", dl.distance); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001); +} + +static void +test_lw_dist2d_arc_arc(void) +{ + /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, + const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, + DISTPTS *dl) */ + 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; + B3.x = 1 ; B3.y = 0; + + /* Arc above the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -1; A1.y = 3; + A2.x = 0 ; A2.y = 2; + A3.x = 1 ; A3.y = 3; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001); + + /* Arc grazes the unit semicircle */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -1; A1.y = 2; + A2.x = 0 ; A2.y = 1; + A3.x = 1 ; A3.y = 2; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001); + + /* Circles intersect, but arcs do not */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -1; A1.y = 1; + A2.x = 0; A2.y = 2; + A3.x = 1; A3.y = 1; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001); + + /* Circles and arcs intersect */ + lw_dist2d_distpts_init(&dl, DIST_MIN); + A1.x = -1; A1.y = 1; + A2.x = 0; A2.y = 0; + A3.x = 1; A3.y = 1; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001); + + /* inscribed and closest on arcs */ + 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; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + //printf("distance %g\n", dl.distance); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001); + + /* inscribed and closest not on arcs */ + 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; + rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl); + //printf("distance %g\n", dl.distance); + CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001); +} + /* ** Used by test harness to register the tests in this file. */ @@ -355,6 +552,9 @@ CU_TestInfo measures_tests[] = PG_TEST(test_rect_tree_intersects_tree), PG_TEST(test_lwgeom_segmentize2d), PG_TEST(test_lwgeom_locate_along), + PG_TEST(test_lw_dist2d_pt_arc), + PG_TEST(test_lw_dist2d_seg_arc), + PG_TEST(test_lw_dist2d_arc_arc), CU_TEST_INFO_NULL }; CU_SuiteInfo measures_suite = {"PostGIS Measures Suite", NULL, NULL, measures_tests}; diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c index 46323432e..72f9456b7 100644 --- a/liblwgeom/cunit/cu_tree.c +++ b/liblwgeom/cunit/cu_tree.c @@ -191,9 +191,9 @@ static void test_tree_circ_distance(void) d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); // printf("d1 = %g d2 = %g\n", d1 * WGS84_RADIUS, d2 * WGS84_RADIUS); // printf("line\n"); - circ_tree_print(c1, 0); +// circ_tree_print(c1, 0); // printf("poly\n"); - circ_tree_print(c2, 0); +// circ_tree_print(c2, 0); circ_tree_free(c1); circ_tree_free(c2); lwgeom_free(lwg1); diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index 10aec67ef..dc32122fd 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -213,7 +213,8 @@ int gbox_overlaps(const GBOX *g1, const GBOX *g2) return LW_TRUE; } -int gbox_overlaps_2d(const GBOX *g1, const GBOX *g2) +int +gbox_overlaps_2d(const GBOX *g1, const GBOX *g2) { /* Make sure our boxes are consistent */ @@ -323,13 +324,13 @@ size_t gbox_serialized_size(uint8_t flags) static int lwcircle_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox) { POINT2D xmin, ymin, xmax, ymax; - POINT4D center; + POINT2D center; int p2_side; double radius; LWDEBUG(2, "lwcircle_calculate_gbox called."); - radius = lwcircle_center(p1, p2, p3, ¢er); + radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, ¢er); /* Negative radius signals straight line, p1/p2/p3 are colinear */ if (radius < 0.0) diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index 9bfc6ceef..b1c7a01ac 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -1802,7 +1802,7 @@ extern char *lwmessage_truncate(char *str, int startpos, int endpos, int maxleng ******************************************************************************/ int lwgeom_has_arc(const LWGEOM *geom); -double lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result); +double lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result); LWGEOM *lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad); LWGEOM *lwgeom_desegmentize(LWGEOM *geom); diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c index 52605f54c..1019d5251 100644 --- a/liblwgeom/lwsegmentize.c +++ b/liblwgeom/lwsegmentize.c @@ -82,13 +82,13 @@ lwgeom_has_arc(const LWGEOM *geom) * point is coincident with either end point, they are taken as colinear. */ double -lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result) +lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result) { - POINT4D c; + POINT2D c; double cx, cy, cr; double temp, bc, cd, det; - c.x = c.y = c.z = c.m = 0.0; + c.x = c.y = 0.0; LWDEBUGF(2, "lwcircle_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y); @@ -157,7 +157,7 @@ static double interpolate_arc(double angle, double a1, double a2, double a3, dou static POINTARRAY * lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad) { - POINT4D center; + POINT2D center; POINT4D pt; int p2_side = 0; int clockwise = LW_TRUE; @@ -170,7 +170,7 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad) LWDEBUG(2, "lwcircle_calculate_gbox called."); - radius = lwcircle_center(p1, p2, p3, ¢er); + radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, ¢er); p2_side = signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2)); /* Matched start/end points imply circle */ @@ -527,16 +527,15 @@ lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad) */ static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b) { - POINT4D center; - POINT4D *centerptr=¢er; - double radius = lwcircle_center(a1, a2, a3, ¢er); + POINT2D center; + double radius = lwcircle_center((POINT2D*)a1, (POINT2D*)a2, (POINT2D*)a3, ¢er); double b_distance, diff; /* Co-linear a1/a2/a3 */ if ( radius < 0.0 ) return LW_FALSE; - b_distance = distance2d_pt_pt((POINT2D*)b, (POINT2D*)centerptr); + b_distance = distance2d_pt_pt((POINT2D*)b, ¢er); diff = fabs(radius - b_distance); LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius); diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index de7765ee4..fee34f628 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -22,11 +22,25 @@ Initializing functions The functions starting the distance-calculation processses --------------------------------------------------------------------------------------------------------------*/ +void +lw_dist2d_distpts_init(DISTPTS *dl, int mode) +{ + dl->twisted = -1; + dl->p1.x = dl->p1.y = 0.0; + dl->p2.x = dl->p2.y = 0.0; + dl->mode = mode; + dl->tolerance = 0.0; + if ( mode == DIST_MIN ) + dl->distance = MAXFLOAT; + else + dl->distance = -1 * MAXFLOAT; +} + /** Function initializing shortestline and longestline calculations. */ LWGEOM * -lw_dist2d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode) +lw_dist2d_distanceline(LWGEOM *lw1, LWGEOM *lw2, int srid, int mode) { double x1,x2,y1,y2; @@ -812,13 +826,417 @@ lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl) } /** +* Returns true if P is on the same side of the plane partition +* defined by A1/A3 as A2 is. Only makes sense if P has already been +* determined to be on the circle defined by A1/A2/A3. +*/ +static int +lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) +{ + return signum(lw_segment_side(A1, A3, A2)) == signum(lw_segment_side(A1, A3, P)); +} + +/** +* Returns true if P is between A1/A2. Only makes sense if P has already been +* deterined to be on the line defined by A1/A2. +*/ +static int +lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2) +{ + return ((A1->x < P->x && P->x < A2->x) || (A1->x > P->x && P->x > A2->x)) || + ((A1->y < P->y && P->y < A2->y) || (A1->y > P->y && P->y > A2->y)); +} + +/** +* Returns true if arc A is actually a point (all vertices are the same) . +*/ +static int +lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) +{ + if ( A1->x == A2->x && A2->x == A3->x && + A1->y == A2->y && A2->y == A3->y ) + return LW_TRUE; + else + return LW_FALSE; +} + +/** +* Calculate the shortest distance between an arc and an edge. +* Line/circle approach from http://stackoverflow.com/questions/1073336/circle-line-collision-detection +*/ +int +lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) +{ + POINT2D C; /* center of arc circle */ + double radius_C; /* radius of arc circle */ + POINT2D D; /* point on A closest to C */ + double dist_C_D; /* distance from C to D */ + int pt_in_arc, pt_in_seg; + DISTPTS dltmp; + + /* Bail out on crazy modes */ + if ( dl->mode < 0 ) + lwerror("lw_dist2d_seg_arc does not support maxdistance mode"); + + /* What if the "arc" is a point? */ + if ( lw_arc_is_pt(B1, B2, B3) ) + return lw_dist2d_pt_seg(B1, A1, A2, dl); + + /* Calculate center and radius of the circle. */ + radius_C = lwcircle_center(B1, B2, B3, &C); + + /* This "arc" is actually a line (B2 is colinear with B1,B3) */ + if ( radius_C < 0.0 ) + return lw_dist2d_seg_seg(A1, A2, B1, B3, dl); + + /* Calculate distance between the line and circle center */ + lw_dist2d_distpts_init(&dltmp, DIST_MIN); + if ( lw_dist2d_pt_seg(&C, A1, A2, &dltmp) == LW_FALSE ) + lwerror("lw_dist2d_pt_seg failed in lw_dist2d_seg_arc"); + + D = dltmp.p1; + dist_C_D = dltmp.distance; + + /* Line intersects circle, maybe arc intersects edge? */ + /* If so, that's the closest point. */ + /* If not, the closest point is one of the end points of A */ + if ( dist_C_D < radius_C ) + { + double length_A; /* length of the segment A */ + POINT2D E, F; /* points of interection of edge A and circle(B) */ + double dist_D_EF; /* distance from D to E or F (same distance both ways) */ + + dist_D_EF = sqrt(radius_C*radius_C - dist_C_D*dist_C_D); + length_A = sqrt((A2->x-A1->x)*(A2->x-A1->x)+(A2->y-A1->y)*(A2->y-A1->y)); + + /* Point of intersection E */ + E.x = D.x - (A2->x-A1->x) * dist_D_EF / length_A; + E.y = D.y - (A2->y-A1->y) * dist_D_EF / length_A; + /* Point of intersection F */ + F.x = D.x + (A2->x-A1->x) * dist_D_EF / length_A; + F.y = D.y + (A2->y-A1->y) * dist_D_EF / length_A; + + + /* If E is within A and within B then it's an interesction point */ + pt_in_arc = lw_pt_in_arc(&E, B1, B2, B3); + pt_in_seg = lw_pt_in_seg(&E, A1, A2); + + if ( pt_in_arc && pt_in_seg ) + { + dl->distance = 0.0; + dl->p1 = E; + dl->p2 = E; + return LW_TRUE; + } + + /* If F is within A and within B then it's an interesction point */ + pt_in_arc = lw_pt_in_arc(&F, B1, B2, B3); + pt_in_seg = lw_pt_in_seg(&F, A1, A2); + + if ( pt_in_arc && pt_in_seg ) + { + dl->distance = 0.0; + dl->p1 = F; + dl->p2 = F; + return LW_TRUE; + } + } + + /* Line grazes circle, maybe arc intersects edge? */ + /* If so, grazing point is the closest point. */ + /* If not, the closest point is one of the end points of A */ + else if ( dist_C_D == radius_C ) + { + /* Closest point D is also the point of grazing */ + pt_in_arc = lw_pt_in_arc(&D, B1, B2, B3); + pt_in_seg = lw_pt_in_seg(&D, A1, A2); + + /* Is D contained in both A and B? */ + if ( pt_in_arc && pt_in_seg ) + { + dl->distance = 0.0; + dl->p1 = D; + dl->p2 = D; + return LW_TRUE; + } + } + /* Line misses circle. */ + /* If closest point to A on circle is within B, then that's the closest */ + /* Otherwise, the closest point will be an end point of A */ + else + { + POINT2D G; /* Point on circle closest to A */ + G.x = C.x + (D.x-C.x) * radius_C / dist_C_D; + G.y = C.y + (D.y-C.y) * radius_C / dist_C_D; + + pt_in_arc = lw_pt_in_arc(&G, B1, B2, B3); + pt_in_seg = lw_pt_in_seg(&D, A1, A2); + + /* Closest point is on the interior of A and B */ + if ( pt_in_arc && pt_in_seg ) + return lw_dist2d_pt_pt(&D, &G, dl); + + } + + /* Now we test the many combinations of end points with either */ + /* arcs or edges. Each previous check determined if the closest */ + /* potential point was within the arc/segment inscribed on the */ + /* line/circle holding the arc/segment. */ + + /* Closest point is in the arc, but not in the segment, so */ + /* one of the segment end points must be the closest. */ + if ( pt_in_arc & ! pt_in_seg ) + { + lw_dist2d_pt_arc(A1, B1, B2, B3, dl); + lw_dist2d_pt_arc(A2, B1, B2, B3, dl); + return LW_TRUE; + } + /* or, one of the arc end points is the closest */ + else if ( pt_in_seg && ! pt_in_arc ) + { + lw_dist2d_pt_seg(B1, A1, A2, dl); + lw_dist2d_pt_seg(B3, A1, A2, dl); + return LW_TRUE; + } + /* Finally, one of the end-point to end-point combos is the closest. */ + else + { + 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); + return LW_TRUE; + } + + return LW_FALSE; +} +int +lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) +{ + double radius_A, d; + POINT2D C; /* center of circle defined by arc A */ + POINT2D X; /* point circle(A) where line from C to P crosses */ + + if ( dl->mode < 0 ) + lwerror("lw_dist2d_pt_arc does not support maxdistance mode"); + + /* What if the arc is a point? */ + if ( lw_arc_is_pt(A1, A2, A3) ) + return lw_dist2d_pt_pt(P, A1, dl); + + /* Calculate centers and radii of circles. */ + radius_A = lwcircle_center(A1, A2, A3, &C); + + /* This "arc" is actually a line (A2 is colinear with A1,A3) */ + if ( radius_A < 0.0 ) + return lw_dist2d_pt_seg(P, A1, A3, dl); + + /* Distance from point to center */ + d = distance2d_pt_pt(&C, P); + + /* X is the point on the circle where the line from P to C crosses */ + X.x = C.x + (P->x - C.x) * radius_A / d; + X.y = C.y + (P->y - C.y) * radius_A / d; + + /* Is crossing point inside the arc? Or arc is actually circle? */ + if ( p2d_same(A1, A3) || lw_pt_in_arc(&X, A1, A2, A3) ) + { + lw_dist2d_pt_pt(P, &X, dl); + } + else + { + /* Distance is the minimum of the distances to the arc end points */ + lw_dist2d_pt_pt(A1, P, dl); + lw_dist2d_pt_pt(A3, P, dl); + } + return LW_TRUE; +} + + +int +lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, + const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, + DISTPTS *dl) +{ + 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 */ + + if ( dl->mode != DIST_MIN ) + lwerror("lw_dist2d_arc_arc only supports mindistance"); + + /* TODO: Handle case where arc is closed circle (A1 = A3) */ + + /* What if one or both of our "arcs" is actually a point? */ + if ( lw_arc_is_pt(B1, B2, B3) && lw_arc_is_pt(A1, A2, A3) ) + return lw_dist2d_pt_pt(B1, A1, dl); + else if ( lw_arc_is_pt(B1, B2, B3) ) + return lw_dist2d_pt_arc(B1, A1, A2, A3, dl); + else if ( lw_arc_is_pt(A1, A2, A3) ) + return lw_dist2d_pt_arc(A1, B1, B2, B3, dl); + + /* Calculate centers and radii of circles. */ + radius_A = lwcircle_center(A1, A2, A3, &CA); + radius_B = lwcircle_center(B1, B2, B3, &CB); + + /* Two co-linear arcs?!? That's two segments. */ + if ( radius_A < 0 && radius_B < 0 ) + return lw_dist2d_seg_seg(A1, A3, B1, B3, dl); + + /* A is co-linear, delegate to lw_dist_seg_arc here. */ + if ( radius_A < 0 ) + return lw_dist2d_seg_arc(A1, A3, B1, B2, B3, dl); + + /* B is co-linear, delegate to lw_dist_seg_arc here. */ + if ( radius_B < 0 ) + return lw_dist2d_seg_arc(B1, B3, A1, A2, A3, dl); + + /* Make sure that arc "A" has the bigger radius */ + if ( radius_B > radius_A ) + { + const POINT2D *tmp; + 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; + } + + /* 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) ) + { + D.x = CA.x + (CB.x - CA.x) * radius_A / d; + D.y = CA.y + (CB.y - CA.y) * radius_A / d; + + pt_in_arc_A = lw_pt_in_arc(&D, A1, A2, A3); + pt_in_arc_B = lw_pt_in_arc(&D, B1, B2, B3); + + /* Arcs do touch at D, return it */ + if ( pt_in_arc_A && pt_in_arc_B ) + { + dl->distance = 0.0; + dl->p1 = D; + dl->p2 = D; + return LW_TRUE; + } + } + /* Disjoint or contained circles don't intersect. Closest point may be on */ + /* the line joining CA to CB. */ + else if ( d > (radius_A + radius_B) /* Disjoint */ || d < (radius_A - radius_B) /* Contained */ ) + { + POINT2D XA, XB; /* Points where the line from CA to CB cross their circle bounds */ + + /* Calculate hypothetical nearest points, the places on the */ + /* two circles where the center-center line crosses. If both */ + /* arcs contain their hypothetical points, that's the crossing distance */ + XA.x = CA.x + (CB.x - CA.x) * radius_A / d; + XA.y = CA.y + (CB.y - CA.y) * radius_A / d; + XB.x = CB.x + (CA.x - CB.x) * radius_B / d; + XB.y = CB.y + (CA.y - CB.y) * radius_B / d; + + pt_in_arc_A = lw_pt_in_arc(&XA, A1, A2, A3); + pt_in_arc_B = lw_pt_in_arc(&XB, B1, B2, B3); + + /* If the nearest points are both within the arcs, that's our answer */ + /* the shortest distance is at the nearest points */ + if ( pt_in_arc_A && pt_in_arc_B ) + { + return lw_dist2d_pt_pt(&XA, &XB, dl); + } + } + /* Circles cross at two points, are either of those points in both arcs? */ + /* http://paulbourke.net/geometry/2circle/ */ + else if ( d < (radius_A + radius_B) ) + { + POINT2D E, F; /* Points where circle(A) and circle(B) cross */ + /* Distance from CA to D */ + double a = (radius_A*radius_A - radius_B*radius_B + d*d) / (2*d); + /* Distance from D to E or F */ + double h = sqrt(radius_A*radius_A - a*a); + + /* Location of D */ + D.x = CA.x + (CB.x - CA.x) * a / d; + D.y = CA.y + (CB.y - CA.y) * a / d; + + /* Start from D and project h units perpendicular to CA-D to get E */ + E.x = D.x + (D.y - CA.y) * h / a; + E.y = D.y + (D.x - CA.x) * h / a; + + /* Crossing point E contained in arcs? */ + pt_in_arc_A = lw_pt_in_arc(&E, A1, A2, A3); + pt_in_arc_B = lw_pt_in_arc(&E, B1, B2, B3); + + if ( pt_in_arc_A && pt_in_arc_B ) + { + dl->p1 = dl->p2 = E; + dl->distance = 0.0; + return LW_TRUE; + } + + /* Start from D and project h units perpendicular to CA-D to get F */ + F.x = D.x - (D.y - CA.y) * h / a; + F.y = D.y - (D.x - CA.x) * h / a; + + /* Crossing point F contained in arcs? */ + pt_in_arc_A = lw_pt_in_arc(&F, A1, A2, A3); + pt_in_arc_B = lw_pt_in_arc(&F, B1, B2, B3); + + if ( pt_in_arc_A && pt_in_arc_B ) + { + dl->p1 = dl->p2 = F; + dl->distance = 0.0; + return LW_TRUE; + } + } + + /* Closest point is in the arc A, but not in the arc B, so */ + /* one of the B end points must be the closest. */ + if ( pt_in_arc_A & ! pt_in_arc_B ) + { + lw_dist2d_pt_arc(B1, A1, A2, A3, dl); + lw_dist2d_pt_arc(B3, A1, A2, A3, dl); + return LW_TRUE; + } + /* Closest point is in the arc B, but not in the arc A, so */ + /* one of the A end points must be the closest. */ + else if ( pt_in_arc_B && ! pt_in_arc_A ) + { + lw_dist2d_pt_arc(A1, B1, B2, B3, dl); + lw_dist2d_pt_arc(A3, B1, B2, B3, dl); + return LW_TRUE; + } + /* Finally, one of the end-point to end-point combos is the closest. */ + else + { + 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); + return LW_TRUE; + } + + return LW_TRUE; +} + +/** Finds the shortest distance between two segments. This function is changed so it is not doing any comparasion of distance but just sending every possible combination further to lw_dist2d_pt_seg */ int -lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl) +lw_dist2d_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl) { double s_top, s_bot,s; double r_top, r_bot,r; @@ -860,8 +1278,8 @@ lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl) If the numerator in eqn 1 is also zero, AB & CD are collinear. */ - r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y) ; - r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x) ; + r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y); + r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x); s_top = (A->y-C->y)*(B->x-A->x) - (A->x-C->x)*(B->y-A->y); s_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x); @@ -1272,7 +1690,7 @@ and just returning the distance without identifying the points. To get this points it was nessecary to change and it also showed to be about 10%faster. */ int -lw_dist2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl) +lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl) { POINT2D c; double r; @@ -1348,7 +1766,7 @@ or most far away from each other depending on dl->mode (max or min) */ int -lw_dist2d_pt_pt(POINT2D *thep1, POINT2D *thep2,DISTPTS *dl) +lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl) { double hside = thep2->x - thep1->x; double vside = thep2->y - thep1->y; @@ -1419,10 +1837,6 @@ distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2) return sqrt ( hside*hside + vside*vside ); - /* the above is more readable - return sqrt( - (p2->x-p1->x) * (p2->x-p1->x) + (p2->y-p1->y) * (p2->y-p1->y) - ); */ } diff --git a/liblwgeom/measures.h b/liblwgeom/measures.h index af6973192..a705f15c9 100644 --- a/liblwgeom/measures.h +++ b/liblwgeom/measures.h @@ -36,18 +36,17 @@ typedef struct /* -Preprocessing functions +* Preprocessing functions */ int lw_dist2d_comp(LWGEOM *lw1, LWGEOM *lw2, DISTPTS *dl); int lw_dist2d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl); int lw_dist2d_check_overlap(LWGEOM *lwg1,LWGEOM *lwg2); int lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); + /* -Brute force functions +* Brute force functions */ - -int lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl); int lw_dist2d_pt_ptarray(POINT2D *p, POINTARRAY *pa, DISTPTS *dl); int lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl); int lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl); @@ -57,16 +56,23 @@ int lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); int lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl); int lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); int lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl); + /* -New faster distance calculations +* New faster distance calculations */ - int lw_dist2d_pre_seg_seg(POINTARRAY *l1, POINTARRAY *l2,LISTSTRUCT *list1, LISTSTRUCT *list2,double k, DISTPTS *dl); int lw_dist2d_selected_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl); int struct_cmp_by_measure(const void *a, const void *b); int lw_dist2d_fast_ptarray_ptarray(POINTARRAY *l1,POINTARRAY *l2, DISTPTS *dl, GBOX *box1, GBOX *box2); + /* -Functions in common for Brute force and new calculation +* Distance calculation primitives. */ -int lw_dist2d_pt_pt(POINT2D *p1, POINT2D *p2, DISTPTS *dl); -int lw_dist2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl); +int lw_dist2d_pt_pt (const POINT2D *P, const POINT2D *Q, DISTPTS *dl); +int lw_dist2d_pt_seg (const POINT2D *P, const POINT2D *A1, const POINT2D *A2, DISTPTS *dl); +int lw_dist2d_pt_arc (const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl); +int lw_dist2d_seg_seg(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, DISTPTS *dl); +int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl); +int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D* B3, DISTPTS *dl); +void lw_dist2d_distpts_init(DISTPTS *dl, int mode); +