From 635f217ae70f0f68014080691bad6e84004124ae Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 23 Oct 2012 22:43:02 +0000 Subject: [PATCH] (#2063) fix the vertex-crossing logic in the circular tree code to use the new edge_intersects routine git-svn-id: http://svn.osgeo.org/postgis/trunk@10532 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/lwgeodetic.c | 13 +------------ liblwgeom/lwgeodetic.h | 14 ++++++++++++++ liblwgeom/lwgeodetic_tree.c | 19 ++++++++++++++----- 3 files changed, 29 insertions(+), 17 deletions(-) diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 9e812b3bc..faf2540a7 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -3123,22 +3123,11 @@ dot_product_side(const POINT3D *p, const POINT3D *q) return dp < 0.0 ? -1 : 1; } -/** -* Bitmask elements for edge_intersects() return value. -*/ -#define PIR_NO_INTERACT 0x00 -#define PIR_INTERSECTS 0x01 -#define PIR_COLINEAR 0x02 -#define PIR_A_TOUCH_LEFT 0x04 -#define PIR_A_TOUCH_RIGHT 0x08 -#define PIR_B_TOUCH_LEFT 0x10 -#define PIR_B_TOUCH_RIGHT 0x20 - /** * Returns non-zero if edges A and B interact. The type of interaction is given in the * return value with the bitmask elements defined above. */ -static int +int edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2) { POINT3D A3, B3; /* Extra vector to build more cross-product-friendly angles */ diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 290b5abaf..0816df712 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -65,6 +65,19 @@ typedef struct */ #define signum(a) ((a) < 0 ? -1 : ((a) > 0 ? 1 : (a))) + +/** +* Bitmask elements for edge_intersects() return value. +*/ +#define PIR_NO_INTERACT 0x00 +#define PIR_INTERSECTS 0x01 +#define PIR_COLINEAR 0x02 +#define PIR_A_TOUCH_LEFT 0x04 +#define PIR_A_TOUCH_RIGHT 0x08 +#define PIR_B_TOUCH_LEFT 0x10 +#define PIR_B_TOUCH_RIGHT 0x20 + + /* * Geodetic calculations */ @@ -86,6 +99,7 @@ int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, G int edge_calculate_gbox(const GEOGRAPHIC_EDGE *e, GBOX *gbox); int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox); int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g); +int edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2); double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest); double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2); void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g); diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c index 94804fcc1..ea57fcfed 100644 --- a/liblwgeom/lwgeodetic_tree.c +++ b/liblwgeom/lwgeodetic_tree.c @@ -394,14 +394,17 @@ circ_nodes_merge(CIRC_NODE** nodes, int num_nodes) */ int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POINT2D* pt_outside, int* on_boundary) { - GEOGRAPHIC_POINT closest, crossing; + GEOGRAPHIC_POINT closest; GEOGRAPHIC_EDGE stab_edge, edge; + POINT3D S1, S2, E1, E2; double d; int i, c; /* Construct a stabline edge from our "inside" to our known outside point */ geographic_point_init(pt->x, pt->y, &(stab_edge.start)); geographic_point_init(pt_outside->x, pt_outside->y, &(stab_edge.end)); + geog2cart(&(stab_edge.start), &S1); + geog2cart(&(stab_edge.end), &S2); LWDEBUG(3, "entered"); @@ -421,18 +424,24 @@ int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POI /* Return the crossing number of this leaf */ if ( circ_node_is_leaf(node) ) { + int inter; LWDEBUGF(3, "leaf node calculation (edge %d)", node->edge_num); geographic_point_init(node->p1->x, node->p1->y, &(edge.start)); geographic_point_init(node->p2->x, node->p2->y, &(edge.end)); - if ( edge_intersection(&stab_edge, &edge, &crossing) ) + geog2cart(&(edge.start), &E1); + geog2cart(&(edge.end), &E2); + + inter = edge_intersects(&S1, &S2, &E1, &E2); + + if ( inter & PIR_INTERSECTS ) { LWDEBUG(3," got stab line edge_intersection with this edge!"); /* To avoid double counting crossings-at-a-vertex, */ /* always ignore crossings at "lower" ends of edges*/ - if ( (FP_EQUALS(crossing.lon, edge.start.lon) && FP_EQUALS(crossing.lat, edge.start.lat) && (edge.start.lat <= edge.end.lat)) || - (FP_EQUALS(crossing.lon, edge.end.lon) && FP_EQUALS(crossing.lat, edge.end.lat) && (edge.end.lat <= edge.start.lat)) ) + + if ( inter & PIR_B_TOUCH_LEFT || inter & PIR_COLINEAR ) { - LWDEBUG(3," rejecting stab line intersection on 'lower' end point vertex"); + LWDEBUG(3," rejecting stab line grazing by left-side edge"); return 0; } else -- 2.40.0