ST_Intersects(geography) returns incorrect result for pure-crossing cases (line cross...

git-svn-id: http://svn.osgeo.org/postgis/trunk@10194 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c
index b27a82a28a2b7c17252e65bd6e2dd5915e4279c8..46323432eef1b5930c4d65ca361e15673cdea431 100644 (file)
--- 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 8d892ceabb37d2dcf31b525e9098eb342c4fe4f6..94804fcc16cff2346a43ac86e8a47aaee4f74aec 100644 (file)
--- 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 )
                        {