Fix a NaN result leaking into the tree building algorithm (optimized 32 bit code...

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

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 0faa8ca2dd4ba0e8f410c66a73b456d9fc6f7220..6e35fcf9144f863936cb51e4a2aed52ff5f45a7c 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -219,6 +219,7 @@ gbox_angular_width(const GBOX* gbox)
                        magnitude = sqrt(pt_n.x*pt_n.x + pt_n.y*pt_n.y);
                        pt_n.x /= magnitude;
                        pt_n.y /= magnitude;
+                       pt_n.z = 0.0;
 
                        dotprod = pt_n.x*pt[j].x + pt_n.y*pt[j].y;
                        angle = acos(dotprod > 1.0 ? 1.0 : dotprod);
@@ -1056,6 +1057,9 @@ int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, G
        n->lat = lat2;
        n->lon = lon2;
 
+       if ( isnan(lat2) || isnan(lon2) )
+               return LW_FAILURE;
+
        return LW_SUCCESS;
 }
 
@@ -2903,4 +2907,4 @@ lwgeom_nudge_geodetic(LWGEOM *geom)
 
        lwerror("unsupported type (%s) passed to lwgeom_nudge_geodetic", lwtype_name(type));
        return rv;
-}
\ No newline at end of file
+}

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index fa7f6db995baa140ed679957d92464e816d2e8e2..0caea337c97029352df1ac82f89270bcc657bf37 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -13,7 +13,10 @@
 #define _LWGEODETIC_H 1
 
 /* For NAN */
+#ifndef _GNU_SOURCE
 #define _GNU_SOURCE
+#endif
+
 #include <math.h>
 
 #ifndef NAN

diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c
index 7b840264e7cb4582c803f2f136b63108ff3ec603..17a009a0850a89cb64586e768a7961919583f860 100644 (file)
--- a/liblwgeom/lwgeodetic_tree.c
+++ b/liblwgeom/lwgeodetic_tree.c
@@ -123,9 +123,7 @@ circ_center_spherical(const GEOGRAPHIC_POINT* c1, const GEOGRAPHIC_POINT* c2, do
                return LW_FAILURE;
        
        /* Center of new circle is projection from start point, using offset distance*/
-       sphere_project(c1, offset, dir, center);
-
-       return LW_SUCCESS;      
+       return sphere_project(c1, offset, dir, center);
 }
 
 /**
@@ -183,7 +181,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
        double offset1, dist, D, r1, ri;
        int i;
 
-       LWDEBUGF(4, "called with %d nodes", num_nodes);
+       LWDEBUGF(3, "called with %d nodes", num_nodes);
 
        /* Can't do anything w/ empty input */
        if ( num_nodes < 1 )
@@ -202,9 +200,9 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
                dist = sphere_distance(&c1, &(c[i]->center));
                ri = c[i]->radius;
 
-               LWDEBUGF(4, "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);
+               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 ( dist < fabs(r1 - ri) )
+               if ( FP_EQUALS(dist, 0) || dist <= fabs(r1 - ri) )
                {
                        /* new contains next */
                        if ( r1 > ri )
@@ -242,7 +240,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
                                new_radius *= 1.1;
                        }
                }
-               LWDEBUGF(4, "new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);      
+               LWDEBUGF(3, "new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);      
        }
        
        node = lwalloc(sizeof(CIRC_NODE));
@@ -294,7 +292,7 @@ circ_tree_new(const POINTARRAY* pa)
        for ( i = 0; i < num_edges; i++ )
        {
                node = circ_node_leaf_new(pa, i);
-               LWDEBUGF(3,"making new leaf node %d", i);
+               LWDEBUGF(4,"making new leaf node %d", i);
                if ( node ) /* Not zero length? */
                        nodes[j++] = node;
        }
@@ -341,7 +339,7 @@ circ_nodes_merge(CIRC_NODE** nodes, int num_nodes)
                        CIRC_NODE **inodes = lwalloc(sizeof(CIRC_NODE*)*node_size);
                        inodes[0] = nodes[node_size*j];
                        inodes[1] = nodes[node_size*j + 1];
-                       LWDEBUGF(3,"merging nodes %d and %d", node_size*j, node_size*j + 1);                    
+                       LWDEBUGF(5,"merging nodes %d and %d", node_size*j, node_size*j + 1);                    
                        nodes[j++] = circ_node_internal_new(inodes, node_size);
                }
                /* Odd number of children, just copy the last node up a level */
@@ -375,7 +373,7 @@ int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POI
        geographic_point_init(pt->x, pt->y, &(stab_edge.start));
        geographic_point_init(pt_outside->x, pt_outside->y, &(stab_edge.end));
        
-       LWDEBUG(4, "entered");
+       LWDEBUG(3, "entered");
        
        /* 
        * If the stabline doesn't cross within the radius of a node, there's no 
@@ -383,9 +381,9 @@ int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POI
        */
 //     circ_tree_print(node, 0);
                
-       LWDEBUGF(4, "working on node %p, edge_num %d, radius %g, center POINT(%g %g)", node, node->edge_num, node->radius, rad2deg(node->center.lon), rad2deg(node->center.lat));
+       LWDEBUGF(3, "working on node %p, edge_num %d, radius %g, center POINT(%g %g)", node, node->edge_num, node->radius, rad2deg(node->center.lon), rad2deg(node->center.lat));
        d = edge_distance_to_point(&stab_edge, &(node->center), &closest);
-       LWDEBUGF(4, "edge_distance_to_point=%g, node_radius=%g", d, node->radius);
+       LWDEBUGF(3, "edge_distance_to_point=%g, node_radius=%g", d, node->radius);
        if ( FP_LTEQ(d, node->radius) )
        {
                LWDEBUGF(3,"entering this branch (%p)", node);
@@ -393,23 +391,23 @@ 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) )
                {
-                       LWDEBUGF(4, "leaf node calculation (edge %d)", node->edge_num);
+                       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) )
                        {
-                               LWDEBUG(4," got stab line edge_intersection with this edge!");
+                               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)) )
                                {
-                                       LWDEBUG(4,"  rejecting stab line intersection on 'lower' end point vertex");
+                                       LWDEBUG(3,"  rejecting stab line intersection on 'lower' end point vertex");
                                        return 0;
                                }
                                else
                                {
-                                       LWDEBUG(4,"  accepting stab line intersection");
+                                       LWDEBUG(3,"  accepting stab line intersection");
                                        return 1;
                                }
                        }

diff --git a/postgis/geography_measurement_trees.c b/postgis/geography_measurement_trees.c
index 892af0460fddce5f5b022c9656608cac659520ea..1b03e83cfd94bd656b55b1f07f4eac2182479ad6 100644 (file)
--- a/postgis/geography_measurement_trees.c
+++ b/postgis/geography_measurement_trees.c
@@ -84,13 +84,17 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
        GEOGRAPHIC_POINT in_gpoint;
        POINT3D in_point3d;
        POINT4D in_point;
-       
+
+       POSTGIS_DEBUGF(3, "tree1_type=%d, lwgeom2->type=%d", tree1_type, lwgeom2->type);
+
        /* If the tree'ed argument is a polygon, do the P-i-P using the tree-based P-i-P */
        if ( tree1_type == POLYGONTYPE || tree1_type == MULTIPOLYGONTYPE )
        {
+               POSTGIS_DEBUG(3, "tree is a polygon, using tree PiP");
                /* Need a gbox to calculate an outside point */
                if ( LW_FAILURE == gserialized_get_gbox_p(g1, &gbox1) )
                {
+                       POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
                        LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
                        lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
                        lwgeom_free(lwgeom1);
@@ -100,6 +104,7 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
                if ( LW_FAILURE == lwgeom_startpoint(lwgeom2, &in_point) )
                {
                        lwerror("CircTreePIP unable to generate start point for lwgeom %p", lwgeom2);
+                       POSTGIS_DEBUG(3, "unable to generate in_point, CircTreePIP returning FALSE");
                        return LW_FALSE;
                }
        
@@ -110,6 +115,7 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
                /* If the candidate isn't in the tree box, it's not in the tree area */
                if ( ! gbox_contains_point3d(&gbox1, &in_point3d) )
                {
+                       POSTGIS_DEBUG(3, "in_point3d is not inside the tree gbox, CircTreePIP returning FALSE");
                        return LW_FALSE;
                }
                /* The candidate point is in the box, so it *might* be inside the tree */
@@ -122,6 +128,7 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
                        /* Calculate a definitive outside point */
                        gbox_pt_outside(&gbox1, &pt2d_outside);
                        /* Test the candidate point for strict containment */
+                       POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
                        return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL);
                }
        }
@@ -131,12 +138,14 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
        {
                int result;
                LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
+               POSTGIS_DEBUG(3, "tree1 not polygonal, but lwgeom2 is, calculating using lwgeom_covers_lwgeom_sphere");
                result = lwgeom_covers_lwgeom_sphere(lwgeom2, lwgeom1);
                lwfree(lwgeom1);
                return result;
        }
        else
        {
+               POSTGIS_DEBUG(3, "neither tree1 nor lwgeom2 polygonal, so CircTreePIP returning FALSE");
                return LW_FALSE;
        }               
 }