From: Paul Ramsey Date: Mon, 24 Mar 2014 15:25:39 +0000 (+0000) Subject: #2634, should fix geography-tree-distance errors in cases X-Git-Tag: 2.2.0rc1~1184 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=001ab8c8f43cb46376fc437b4f09f312f6bd2fd5;p=postgis #2634, should fix geography-tree-distance errors in cases where polygons interact with collections. git-svn-id: http://svn.osgeo.org/postgis/trunk@12340 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c index a224676d2..92610fd6f 100644 --- a/liblwgeom/cunit/cu_tree.c +++ b/liblwgeom/cunit/cu_tree.c @@ -16,8 +16,8 @@ #include "CUnit/Basic.h" #include "liblwgeom_internal.h" -#include "lwgeodetic.h"; -#include "lwgeodetic_tree.h"; +#include "lwgeodetic.h" +#include "lwgeodetic_tree.h" #include "cu_tester.h" @@ -177,11 +177,12 @@ static void test_tree_circ_distance(void) LWGEOM *lwg1, *lwg2; CIRC_NODE *c1, *c2; SPHEROID s; - double d1, d2; + double d1, d2, d3, d4; double threshold = 0.0; spheroid_init(&s, 1.0, 1.0); +#if 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); @@ -283,6 +284,52 @@ static void test_tree_circ_distance(void) lwgeom_free(lwg1); lwgeom_free(lwg2); CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.0000001); + +#endif + + /* Ticket #2634 */ + lwg1 = lwgeom_from_wkt("MULTIPOINT (-10 40,-10 65,10 40,10 65,30 40,30 65,50 40,50 65)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("POLYGON((-9.1111111 40,-9.14954053919354 39.6098193559677,-9.26335203497743 39.2346331352698,-9.44817187539491 38.8888595339608,-9.6968975376269 38.5857864376269,-9.99997063396079 38.3370607753949,-10.3457442352698 38.1522409349774,-10.7209304559677 38.0384294391935,-11.1111111 38,-11.5012917440323 38.0384294391935,-11.8764779647302 38.1522409349774,-12.2222515660392 38.3370607753949,-12.5253246623731 38.5857864376269,-12.7740503246051 38.8888595339608,-12.9588701650226 39.2346331352698,-13.0726816608065 39.6098193559677,-13.1111111 40,-13.0726816608065 40.3901806440322,-12.9588701650226 40.7653668647302,-12.7740503246051 41.1111404660392,-12.5253246623731 41.4142135623731,-12.2222515660392 41.6629392246051,-11.8764779647302 41.8477590650226,-11.5012917440323 41.9615705608065,-11.1111111 42,-10.7209304559678 41.9615705608065,-10.3457442352698 41.8477590650226,-9.9999706339608 41.6629392246051,-9.69689753762691 41.4142135623731,-9.44817187539491 41.1111404660392,-9.26335203497743 40.7653668647302,-9.14954053919354 40.3901806440323,-9.1111111 40))", 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("multipoint\n"); + circ_tree_print(c1, 0); + // printf("polygon\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); + +#if 0 + /* Ticket #2634 */ + lwg1 = lwgeom_from_wkt("MULTIPOINT Z (-10 40 1,-10 65 1,10 40 1,10 65 1,30 40 1,30 65 1,50 40 1,50 65 1,-10 40 2,-10 65 2,10 40 2,10 65 2,30 40 2,30 65 2,50 40 2,50 65 2,-10 40 3,-10 65 3,10 40 3,10 65 3,30 40 3,30 65 3,50 40 3,50 65 3)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_from_wkt("MULTIPOLYGON(((-9.1111111 40,-9.14954053919354 39.6098193559677,-9.26335203497743 39.2346331352698,-9.44817187539491 38.8888595339608,-9.6968975376269 38.5857864376269,-9.99997063396079 38.3370607753949,-10.3457442352698 38.1522409349774,-10.7209304559677 38.0384294391935,-11.1111111 38,-11.5012917440323 38.0384294391935,-11.8764779647302 38.1522409349774,-12.2222515660392 38.3370607753949,-12.5253246623731 38.5857864376269,-12.7740503246051 38.8888595339608,-12.9588701650226 39.2346331352698,-13.0726816608065 39.6098193559677,-13.1111111 40,-13.0726816608065 40.3901806440322,-12.9588701650226 40.7653668647302,-12.7740503246051 41.1111404660392,-12.5253246623731 41.4142135623731,-12.2222515660392 41.6629392246051,-11.8764779647302 41.8477590650226,-11.5012917440323 41.9615705608065,-11.1111111 42,-10.7209304559678 41.9615705608065,-10.3457442352698 41.8477590650226,-9.9999706339608 41.6629392246051,-9.69689753762691 41.4142135623731,-9.44817187539491 41.1111404660392,-9.26335203497743 40.7653668647302,-9.14954053919354 40.3901806440323,-9.1111111 40)))", LW_PARSER_CHECK_NONE); + c1 = lwgeom_calculate_circ_tree(lwg1); + c2 = lwgeom_calculate_circ_tree(lwg2); + printf("\n"); + circ_tree_print(c1, 0); + printf("\n"); + circ_tree_print(c2, 0); + printf("\n"); + d1 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold); + d2 = circ_tree_distance_tree(c1, c2, &s, threshold); + d3 = circ_tree_distance_tree(c1, c2, &s, threshold); + d4 = circ_tree_distance_tree(c1, c2, &s, threshold); + printf("\n d1-no-tree %20.20g\n d2 %20.20g\n d3 %20.20g\n d4 %20.20g\n", d1, d2, d3, d4); + circ_tree_free(c1); + circ_tree_free(c2); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + CU_ASSERT_DOUBLE_EQUAL(d1, d2, 0.00000001); + CU_ASSERT_DOUBLE_EQUAL(d1, d3, 0.00000001); + CU_ASSERT_DOUBLE_EQUAL(d1, d4, 0.00000001); +#endif + } diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index be5a9255b..93cce9ab0 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1419,7 +1419,22 @@ int edge_calculate_gbox(const POINT3D *A1, const POINT3D *A2, GBOX *gbox) return LW_SUCCESS; } - +void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside) +{ + /* Make sure we have boxes */ + if ( poly->bbox ) + { + gbox_pt_outside(poly->bbox, pt_outside); + return; + } + else + { + GBOX gbox; + lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox); + gbox_pt_outside(&gbox, pt_outside); + return; + } +} /** * Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index acd10251f..74c976d01 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -105,6 +105,7 @@ double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g); int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test); int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test); +void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside); int ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test); double ptarray_area_sphere(const POINTARRAY *pa); double latitude_degrees_normalize(double lat); @@ -130,4 +131,16 @@ double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, c double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid); int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g); + #endif /* _LWGEODETIC_H */ + + + +/** +* Notes for rewrite +* +* Define separate POINT types for 2-d-points-in-radiands and 3-d-points-in-geocentric +* Maintain consistent units (radians?) throughout all calculations +* Put an index pointer onto LWGEOM itself, and cache the indexed LWGEOM instead of a bare tree +* only primitive objects should get a tree +*/ \ No newline at end of file diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c index 383f6ac75..651e41bd7 100644 --- a/liblwgeom/lwgeodetic_tree.c +++ b/liblwgeom/lwgeodetic_tree.c @@ -82,7 +82,8 @@ circ_node_leaf_new(const POINTARRAY* pa, int i) node->edge_num = i; /* Zero out metadata */ - node->pt_outside = NULL; + node->pt_outside.x = 0.0; + node->pt_outside.y = 0.0; node->geom_type = 0; return node; @@ -102,7 +103,8 @@ circ_node_leaf_point_new(const POINTARRAY* pa) tree->num_nodes = 0; tree->edge_num = 0; tree->geom_type = POINTTYPE; - tree->pt_outside = NULL; + tree->pt_outside.x = 0.0; + tree->pt_outside.y = 0.0; return tree; } @@ -205,7 +207,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes) GEOGRAPHIC_POINT new_center, c1; double new_radius; double offset1, dist, D, r1, ri; - int i; + int i, new_geom_type; LWDEBUGF(3, "called with %d nodes --", num_nodes); @@ -216,6 +218,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes) /* Initialize calculation with values of the first circle */ new_center = c[0]->center; new_radius = c[0]->radius; + new_geom_type = c[0]->geom_type; /* Merge each remaining circle into the new circle */ for ( i = 1; i < num_nodes; i++ ) @@ -226,6 +229,32 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes) dist = sphere_distance(&c1, &(c[i]->center)); ri = c[i]->radius; + /* Promote geometry types up the tree, getting more and more collected */ + /* Go until we find a value */ + if ( ! new_geom_type ) + { + new_geom_type = c[i]->geom_type; + } + /* Promote singleton to a multi-type */ + else if ( ! lwtype_is_collection(new_geom_type) ) + { + /* Anonymous collection if types differ */ + if ( new_geom_type != c[i]->geom_type ) + { + new_geom_type = COLLECTIONTYPE; + } + else + { + new_geom_type = lwtype_get_collectiontype(new_geom_type); + } + } + /* If we can't add next feature to this collection cleanly, promote again to anonymous collection */ + else if ( new_geom_type != lwtype_get_collectiontype(c[i]->geom_type) ) + { + new_geom_type = COLLECTIONTYPE; + } + + LWDEBUGF(3, "distance between new (%g %g) and %i (%g %g) is %g", c1.lon, c1.lat, i, c[i]->center.lon, c[i]->center.lat, dist); if ( FP_EQUALS(dist, 0) ) @@ -286,8 +315,9 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes) node->num_nodes = num_nodes; node->nodes = c; node->edge_num = -1; - node->geom_type = 0; - node->pt_outside = NULL; + node->geom_type = new_geom_type; + node->pt_outside.x = 0.0; + node->pt_outside.y = 0.0; return node; } @@ -358,6 +388,9 @@ circ_nodes_merge(CIRC_NODE** nodes, int num_nodes) int num_parents = 0; int j; + // TODO, roll geom_type *up* as tree is built, changing to collection types as simple types are merged + // TODO, change the distance algorithm to drive down to simple types first, test pip on poly/other cases, then test edges + while( num_children > 1 ) { for ( j = 0; j < num_children; j++ ) @@ -397,7 +430,7 @@ circ_nodes_merge(CIRC_NODE** nodes, int num_nodes) /** * Returns a #POINT2D that is a vertex of the input shape */ -int circ_tree_get_point(const CIRC_NODE* node, POINT4D* pt) +int circ_tree_get_point(const CIRC_NODE* node, POINT2D* pt) { if ( circ_node_is_leaf(node) ) { @@ -544,12 +577,12 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl double d, d_min; int i; - LWDEBUGF(4, "entered, min_dist %.8g max_dist %.8g", *min_dist, *max_dist); + LWDEBUGF(4, "entered, min_dist=%.8g max_dist=%.8g, type1=%d, type2=%d", *min_dist, *max_dist, n1->geom_type, n2->geom_type); // circ_tree_print(n1, 0); // circ_tree_print(n2, 0); /* Short circuit if we've already hit the minimum */ - if( FP_LT(*min_dist, threshold) ) + if( *min_dist <= threshold ) return *min_dist; /* If your minimum is greater than anyone's maximum, you can't hold the winner */ @@ -565,127 +598,204 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl if( max < *max_dist ) *max_dist = max; + /* Polygon on one side, primitive type on the other. Check for point-in-polygon */ + /* short circuit. */ + if ( n1->geom_type == POLYGONTYPE && n2->geom_type && ! lwtype_is_collection(n2->geom_type) ) + { + POINT2D pt; + circ_tree_get_point(n2, &pt); + LWDEBUGF(4, "n1 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y); + if ( circ_tree_contains_point(n1, &pt, &(n1->pt_outside), NULL) ) + { + LWDEBUG(4, "it does"); + *min_dist = 0.0; + geographic_point_init(pt.x, pt.y, closest1); + geographic_point_init(pt.x, pt.y, closest2); + return *min_dist; + } + } + /* Polygon on one side, primitive type on the other. Check for point-in-polygon */ + /* short circuit. */ + if ( n2->geom_type == POLYGONTYPE && n1->geom_type && ! lwtype_is_collection(n1->geom_type) ) + { + POINT2D pt; + circ_tree_get_point(n1, &pt); + LWDEBUGF(4, "n2 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y); + if ( circ_tree_contains_point(n2, &pt, &(n2->pt_outside), NULL) ) + { + LWDEBUG(4, "it does"); + geographic_point_init(pt.x, pt.y, closest1); + geographic_point_init(pt.x, pt.y, closest2); + *min_dist = 0.0; + return *min_dist; + } + } + /* Both leaf nodes, do a real distance calculation */ - if( circ_node_is_leaf(n1) ) + if( circ_node_is_leaf(n1) && circ_node_is_leaf(n2) ) { - if( circ_node_is_leaf(n2) ) + double d; + GEOGRAPHIC_POINT close1, close2; + 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 ) { - double d; - GEOGRAPHIC_POINT close1, close2; - 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 ) - { - GEOGRAPHIC_EDGE e; - GEOGRAPHIC_POINT gp1, gp2; + GEOGRAPHIC_EDGE e; + GEOGRAPHIC_POINT gp1, gp2; - /* Both nodes are points! */ - if ( n1->p1 == n1->p2 && n2->p1 == n2->p2 ) - { - geographic_point_init(n1->p1->x, n1->p1->y, &gp1); - geographic_point_init(n2->p1->x, n2->p1->y, &gp2); - close1 = gp1; close2 = gp2; - d = sphere_distance(&gp1, &gp2); - } - /* Node 1 is a point */ - else if ( n1->p1 == n1->p2 ) - { - geographic_point_init(n1->p1->x, n1->p1->y, &gp1); - geographic_point_init(n2->p1->x, n2->p1->y, &(e.start)); - geographic_point_init(n2->p2->x, n2->p2->y, &(e.end)); - close1 = gp1; - d = edge_distance_to_point(&e, &gp1, &close2); - } - /* Node 2 is a point */ - else - { - geographic_point_init(n2->p1->x, n2->p1->y, &gp1); - geographic_point_init(n1->p1->x, n1->p1->y, &(e.start)); - geographic_point_init(n1->p2->x, n1->p2->y, &(e.end)); - close1 = gp1; - d = edge_distance_to_point(&e, &gp1, &close2); - } - LWDEBUGF(4, " got distance %g", d); + /* Both nodes are points! */ + if ( n1->p1 == n1->p2 && n2->p1 == n2->p2 ) + { + geographic_point_init(n1->p1->x, n1->p1->y, &gp1); + geographic_point_init(n2->p1->x, n2->p1->y, &gp2); + close1 = gp1; close2 = gp2; + d = sphere_distance(&gp1, &gp2); + } + /* Node 1 is a point */ + else if ( n1->p1 == n1->p2 ) + { + geographic_point_init(n1->p1->x, n1->p1->y, &gp1); + geographic_point_init(n2->p1->x, n2->p1->y, &(e.start)); + geographic_point_init(n2->p2->x, n2->p2->y, &(e.end)); + close1 = gp1; + d = edge_distance_to_point(&e, &gp1, &close2); } - /* Both nodes are edges */ + /* Node 2 is a point */ else { - GEOGRAPHIC_EDGE e1, e2; - GEOGRAPHIC_POINT g; - POINT3D A1, A2, B1, B2; - 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)); - geog2cart(&(e1.start), &A1); - geog2cart(&(e1.end), &A2); - geog2cart(&(e2.start), &B1); - geog2cart(&(e2.end), &B2); - if ( edge_intersects(&A1, &A2, &B1, &B2) ) - { - d = 0.0; - edge_intersection(&e1, &e2, &g); - close1 = close2 = g; - } - else - { - d = edge_distance_to_edge(&e1, &e2, &close1, &close2); - } - LWDEBUGF(4, "edge_distance_to_edge returned %g", d); + geographic_point_init(n2->p1->x, n2->p1->y, &gp1); + geographic_point_init(n1->p1->x, n1->p1->y, &(e.start)); + geographic_point_init(n1->p2->x, n1->p2->y, &(e.end)); + close1 = gp1; + d = edge_distance_to_point(&e, &gp1, &close2); } - if ( d < *min_dist ) + LWDEBUGF(4, " got distance %g", d); + } + /* Both nodes are edges */ + else + { + GEOGRAPHIC_EDGE e1, e2; + GEOGRAPHIC_POINT g; + POINT3D A1, A2, B1, B2; + 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)); + geog2cart(&(e1.start), &A1); + geog2cart(&(e1.end), &A2); + geog2cart(&(e2.start), &B1); + geog2cart(&(e2.end), &B2); + if ( edge_intersects(&A1, &A2, &B1, &B2) ) { - *min_dist = d; - *closest1 = close1; - *closest2 = close2; + d = 0.0; + edge_intersection(&e1, &e2, &g); + close1 = close2 = g; } - return d; + else + { + d = edge_distance_to_edge(&e1, &e2, &close1, &close2); + } + LWDEBUGF(4, "edge_distance_to_edge returned %g", d); } - else + if ( d < *min_dist ) + { + *min_dist = d; + *closest1 = close1; + *closest2 = close2; + } + return d; + } + else + { + d_min = MAXFLOAT; + /* Drive the recursion into the COLLECTION types first so we end up with */ + /* pairings of primitive geometries that can be forced into the point-in-polygon */ + /* tests above. */ + if ( n1->geom_type && lwtype_is_collection(n1->geom_type) ) + { + for ( i = 0; i < n1->num_nodes; i++ ) + { + d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2); + d_min = FP_MIN(d_min, d); + } + } + else if ( n2->geom_type && lwtype_is_collection(n2->geom_type) ) { - d_min = MAXFLOAT; for ( i = 0; i < n2->num_nodes; i++ ) { d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2); d_min = FP_MIN(d_min, d); } - return d_min; } - } - else - { - d_min = MAXFLOAT; - for ( i = 0; i < n1->num_nodes; i++ ) + else if ( ! circ_node_is_leaf(n1) ) + { + for ( i = 0; i < n1->num_nodes; i++ ) + { + d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2); + d_min = FP_MIN(d_min, d); + } + } + else if ( ! circ_node_is_leaf(n2) ) { - d = circ_tree_distance_tree_internal(n2, n1->nodes[i], threshold, min_dist, max_dist, closest1, closest2); - d_min = FP_MIN(d_min, d); + for ( i = 0; i < n2->num_nodes; i++ ) + { + d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2); + d_min = FP_MIN(d_min, d); + } } + else + { + /* Never get here */ + } + return d_min; } } + + + void circ_tree_print(const CIRC_NODE* node, int depth) { int i; - - if ( node->num_nodes > 0 ) - { - printf("%*s C(%.5g %.5g) R(%.5g)\n", - 3*depth + 6, "NODE", - node->center.lon, node->center.lat, - node->radius - ); - } - else + + if (circ_node_is_leaf(node)) { - printf("%*s[%d] C(%.5g %.5g) R(%.5g) ((%.5g %.5g),(%.5g,%.5g))\n", + printf("%*s[%d] C(%.5g %.5g) R(%.5g) ((%.5g %.5g),(%.5g,%.5g))", 3*depth + 6, "NODE", node->edge_num, node->center.lon, node->center.lat, node->radius, node->p1->x, node->p1->y, node->p2->x, node->p2->y ); + if ( node->geom_type ) + { + printf(" %s", lwtype_name(node->geom_type)); + } + if ( node->geom_type == POLYGONTYPE ) + { + printf(" O(%.5g %.5g)", node->pt_outside.x, node->pt_outside.y); + } + printf("\n"); + + } + else + { + printf("%*s C(%.5g %.5g) R(%.5g)", + 3*depth + 6, "NODE", + node->center.lon, node->center.lat, + node->radius + ); + if ( node->geom_type ) + { + printf(" %s", lwtype_name(node->geom_type)); + } + if ( node->geom_type == POLYGONTYPE ) + { + printf(" O(%.5g %.5g)", node->pt_outside.x, node->pt_outside.y); + } + printf("\n"); } for ( i = 0; i < node->num_nodes; i++ ) { @@ -722,24 +832,31 @@ lwpoly_calculate_circ_tree(const LWPOLY* lwpoly) /* One ring? Handle it like a line. */ if ( lwpoly->nrings == 1 ) - return circ_tree_new(lwpoly->rings[0]); - - /* Calculate a tree for each non-trivial ring of the polygon */ - nodes = lwalloc(lwpoly->nrings * sizeof(CIRC_NODE*)); - for ( i = 0; i < lwpoly->nrings; i++ ) { - node = circ_tree_new(lwpoly->rings[i]); - if ( node ) - nodes[j++] = node; + node = circ_tree_new(lwpoly->rings[0]); } - /* Put the trees into a spatially correlated order */ - circ_nodes_sort(nodes, j); - /* Merge the trees pairwise up to a parent node and return */ - node = circ_nodes_merge(nodes, j); - /* Don't need the working list any more */ - lwfree(nodes); - + else + { + /* Calculate a tree for each non-trivial ring of the polygon */ + nodes = lwalloc(lwpoly->nrings * sizeof(CIRC_NODE*)); + for ( i = 0; i < lwpoly->nrings; i++ ) + { + node = circ_tree_new(lwpoly->rings[i]); + if ( node ) + nodes[j++] = node; + } + /* Put the trees into a spatially correlated order */ + circ_nodes_sort(nodes, j); + /* Merge the trees pairwise up to a parent node and return */ + node = circ_nodes_merge(nodes, j); + /* Don't need the working list any more */ + lwfree(nodes); + } + + /* Metatdata about polygons, we need this to apply P-i-P tests */ + /* selectively when doing distance calculations */ node->geom_type = lwgeom_get_type((LWGEOM*)lwpoly); + lwpoly_pt_outside(lwpoly, &(node->pt_outside)); return node; } @@ -770,7 +887,7 @@ lwcollection_calculate_circ_tree(const LWCOLLECTION* lwcol) /* Don't need the working list any more */ lwfree(nodes); - node->geom_type = lwgeom_get_type((LWGEOM*)lwcol);; + node->geom_type = lwgeom_get_type((LWGEOM*)lwcol); return node; } diff --git a/liblwgeom/lwgeodetic_tree.h b/liblwgeom/lwgeodetic_tree.h index e2878f584..af772bb70 100644 --- a/liblwgeom/lwgeodetic_tree.h +++ b/liblwgeom/lwgeodetic_tree.h @@ -17,7 +17,7 @@ typedef struct circ_node struct circ_node** nodes; int edge_num; int geom_type; - POINT2D* pt_outside; + POINT2D pt_outside; POINT2D* p1; POINT2D* p2; } CIRC_NODE; @@ -28,6 +28,8 @@ void circ_tree_free(CIRC_NODE* node); int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POINT2D* pt_outside, int* on_boundary); double circ_tree_distance_tree(const CIRC_NODE* n1, const CIRC_NODE* n2, const SPHEROID *spheroid, double threshold); CIRC_NODE* lwgeom_calculate_circ_tree(const LWGEOM* lwgeom); -int circ_tree_get_point(const CIRC_NODE* node, POINT4D* pt); +int circ_tree_get_point(const CIRC_NODE* node, POINT2D* pt); #endif /* _LWGEODETIC_TREE_H */ + + diff --git a/postgis/geography_measurement_trees.c b/postgis/geography_measurement_trees.c index 387661abb..59f537ef2 100644 --- a/postgis/geography_measurement_trees.c +++ b/postgis/geography_measurement_trees.c @@ -199,7 +199,10 @@ geography_distance_cache_tolerance(FunctionCallInfoData* fcinfo, const GSERIALIZ circtree = lwgeom_calculate_circ_tree(lwgeom); if ( geomtype == POLYGONTYPE || geomtype == MULTIPOLYGONTYPE ) { - circ_tree_get_point(circtree_cached, &p4d); + POINT2D p2d; + circ_tree_get_point(circtree_cached, &p2d); + p4d.x = p2d.x; + p4d.y = p2d.y; if ( CircTreePIP(circtree, g, &p4d) ) { *distance = 0.0;