From: Paul Ramsey Date: Tue, 10 Jul 2012 20:32:54 +0000 (+0000) Subject: Fix a NaN result leaking into the tree building algorithm (optimized 32 bit code... X-Git-Tag: 2.1.0beta2~814 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=b28a605814eb31ca1f16a3b7818d604a8bdacd88;p=postgis Fix a NaN result leaking into the tree building algorithm (optimized 32 bit code only!). 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 0faa8ca2d..6e35fcf91 100644 --- 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 fa7f6db99..0caea337c 100644 --- 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 #ifndef NAN diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c index 7b840264e..17a009a08 100644 --- 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 892af0460..1b03e83cf 100644 --- 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; } }