(#2063) fix the vertex-crossing logic in the circular tree code to use the new edge_i...

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

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 9e812b3bcd09240a1179e27a29856fa6ff94a6a0..faf2540a749d2ec93b2888949a53eebb2457d51d 100644 (file)
--- 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 290b5abaf0e24834d352ce184d244f174ebf37d7..0816df712538edcdded8502a1113b6f25bda8a73 100644 (file)
--- 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 94804fcc16cff2346a43ac86e8a47aaee4f74aec..ea57fcfed876410609e95e878c01403af3077ec1 100644 (file)
--- 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