#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"
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);
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
+
}
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;
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;
}
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);
/* 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++ )
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) )
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;
}
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++ )
/**
* 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) )
{
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 */
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++ )
{
/* 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;
}
/* 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;
}
projects / postgis / commitdiff