From da9ecad04310d76d74a54467fa61774b814e9829 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 21 Aug 2012 18:22:41 +0000 Subject: [PATCH] ST_Intersects(geography) returns incorrect result for pure-crossing cases (line cross line, line crosses polygon) (#1958) git-svn-id: http://svn.osgeo.org/postgis/trunk@10194 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_tree.c | 143 ++++++++++++++++++++---------------- liblwgeom/lwgeodetic_tree.c | 14 +++- 2 files changed, 92 insertions(+), 65 deletions(-) diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c index b27a82a28..46323432e 100644 --- a/liblwgeom/cunit/cu_tree.c +++ b/liblwgeom/cunit/cu_tree.c @@ -174,79 +174,96 @@ static void test_tree_circ_pip2(void) static void test_tree_circ_distance(void) { - LWLINE *line; - LWPOINT *point; - CIRC_NODE *cline, *cpoint; + LWGEOM *lwg1, *lwg2; + CIRC_NODE *c1, *c2; SPHEROID s; - double distance_tree, distance_geom; + double d1, d2; double threshold = 0.0; spheroid_init(&s, 1.0, 1.0); + /* Ticket #1958 */ + lwg1 = lwgeom_from_wkt("LINESTRING(22.88333 41.96667,21.32667 42.13667)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POLYGON((22.94472 41.34667,22.87528 41.99028,22.87389 41.98472,22.87472 41.98333,22.94472 41.34667))", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + 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); +// printf("poly\n"); + circ_tree_print(c2, 0); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.0000001); + /* Ticket #1951 */ - line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING(0 0, 0 0)", LW_PARSER_CHECK_NONE)); - point = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT(0.1 0.1)", LW_PARSER_CHECK_NONE)); - cline = circ_tree_new(line->points); - cpoint = circ_tree_new(point->point); - distance_tree = circ_tree_distance_tree(cpoint, cline, &s, threshold); - distance_geom = lwgeom_distance_spheroid((LWGEOM*)line, (LWGEOM*)point, &s, threshold); - circ_tree_free(cline); - circ_tree_free(cpoint); - lwline_free(line); - lwpoint_free(point); - CU_ASSERT_DOUBLE_EQUAL(distance_geom, distance_geom, 0.0001); - - line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE)); - point = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT(-2 0)", LW_PARSER_CHECK_NONE)); - cline = circ_tree_new(line->points); - cpoint = circ_tree_new(point->point); - distance_tree = circ_tree_distance_tree(cpoint, cline, &s, threshold); - distance_geom = lwgeom_distance_spheroid((LWGEOM*)line, (LWGEOM*)point, &s, threshold); - circ_tree_free(cline); - circ_tree_free(cpoint); - lwline_free(line); - lwpoint_free(point); - CU_ASSERT_DOUBLE_EQUAL(distance_geom, distance_geom, 0.0001); - - line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE)); - point = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT(2 2)", LW_PARSER_CHECK_NONE)); - cline = circ_tree_new(line->points); - cpoint = circ_tree_new(point->point); - distance_tree = circ_tree_distance_tree(cpoint, cline, &s, threshold); - distance_geom = lwgeom_distance_spheroid((LWGEOM*)line, (LWGEOM*)point, &s, threshold); - circ_tree_free(cline); - circ_tree_free(cpoint); - lwline_free(line); - lwpoint_free(point); - CU_ASSERT_DOUBLE_EQUAL(distance_geom, distance_geom, 0.0001); - - line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE)); - point = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT(1 1)", LW_PARSER_CHECK_NONE)); - cline = circ_tree_new(line->points); - cpoint = circ_tree_new(point->point); - distance_tree = circ_tree_distance_tree(cpoint, cline, &s, threshold); - distance_geom = lwgeom_distance_spheroid((LWGEOM*)line, (LWGEOM*)point, &s, threshold); - circ_tree_free(cline); - circ_tree_free(cpoint); - lwline_free(line); - lwpoint_free(point); - CU_ASSERT_DOUBLE_EQUAL(distance_tree, distance_geom, 0.0001); - - line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE)); - point = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT(1 0.5)", LW_PARSER_CHECK_NONE)); - cline = circ_tree_new(line->points); - cpoint = circ_tree_new(point->point); - distance_tree = circ_tree_distance_tree(cpoint, cline, &s, threshold); - distance_geom = lwgeom_distance_spheroid((LWGEOM*)line, (LWGEOM*)point, &s, threshold); + lwg1 = lwgeom_from_wkt("LINESTRING(0 0, 0 0)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POINT(0.1 0.1)", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00001); + + lwg1 = lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POINT(-2 0)", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00001); + + lwg1 = lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POINT(2 2)", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00001); + + lwg1 = lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POINT(1 1)", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00001); + + lwg1 = lwgeom_from_wkt("LINESTRING(-1 -1,0 -1,1 -1,1 0,1 1,0 0,-1 1,-1 0,-1 -1)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POINT(1 0.5)", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + d1 = circ_tree_distance_tree(c1, c2, &s, threshold); + d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); // printf("distance_tree %g\n", distance_tree); // printf("distance_geom %g\n", distance_geom); // circ_tree_print(cline, 0); // circ_tree_print(cpoint, 0); - circ_tree_free(cline); - circ_tree_free(cpoint); - lwline_free(line); - lwpoint_free(point); - CU_ASSERT_DOUBLE_EQUAL(distance_tree, distance_geom, 0.0001); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00001); } diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c index 8d892ceab..94804fcc1 100644 --- a/liblwgeom/lwgeodetic_tree.c +++ b/liblwgeom/lwgeodetic_tree.c @@ -532,7 +532,7 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl { double d; GEOGRAPHIC_POINT close1, close2; - LWDEBUGF(4, "testing pair [%d], [%d]", n1->edge_num, n2->edge_num); + LWDEBUGF(4, "testing leaf pair [%d], [%d]", n1->edge_num, n2->edge_num); /* One of the nodes is a point */ if ( n1->p1 == n1->p2 || n2->p1 == n2->p2 ) { @@ -571,11 +571,21 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl else { GEOGRAPHIC_EDGE e1, e2; + GEOGRAPHIC_POINT g1; geographic_point_init(n1->p1->x, n1->p1->y, &(e1.start)); geographic_point_init(n1->p2->x, n1->p2->y, &(e1.end)); geographic_point_init(n2->p1->x, n2->p1->y, &(e2.start)); geographic_point_init(n2->p2->x, n2->p2->y, &(e2.end)); - d = edge_distance_to_edge(&e1, &e2, &close1, &close2); + if ( edge_intersection(&e1, &e2, &g1) ) + { + d = 0.0; + close1 = close2 = g1; + } + else + { + d = edge_distance_to_edge(&e1, &e2, &close1, &close2); + } + LWDEBUGF(4, "edge_distance_to_edge returned %g", d); } if ( d < *min_dist ) { -- 2.40.0