From ab66d2ddf1ee79eb2c04e6b448a97153ce1fad85 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 4 Feb 2014 19:44:00 +0000 Subject: [PATCH] #2556, geography ST_Intersects results depending on insert order git-svn-id: http://svn.osgeo.org/postgis/trunk@12219 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + postgis/geography_measurement_trees.c | 190 +++++++++++--------------- regress/tickets.sql | 9 ++ regress/tickets_expected | 2 + 4 files changed, 91 insertions(+), 111 deletions(-) diff --git a/NEWS b/NEWS index 8136d9d61..98f617560 100644 --- a/NEWS +++ b/NEWS @@ -48,6 +48,7 @@ PostGIS 2.2.0 - #2555, Fix parsing issue of range arguments of ST_Reclass - #2589, Remove use of unnecessary void pointers - #2607, Cannot open more than 1024 out-db files in one process + - #2556, geography ST_Intersects results depending on insert order * Code refactoring * diff --git a/postgis/geography_measurement_trees.c b/postgis/geography_measurement_trees.c index ed163ea9c..a80ab5da8 100644 --- a/postgis/geography_measurement_trees.c +++ b/postgis/geography_measurement_trees.c @@ -9,7 +9,7 @@ * liblwgeom, where most of the circtree logic lives. */ typedef struct { - int type; // + int type; // GSERIALIZED* geom1; // GSERIALIZED* geom2; // size_t geom1_size; // @@ -76,14 +76,14 @@ GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2); } + static int -CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2) +CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point) { int tree1_type = gserialized_get_type(g1); GBOX gbox1; GEOGRAPHIC_POINT in_gpoint; POINT3D in_point3d; - POINT4D in_point; POSTGIS_DEBUGF(3, "tree1_type=%d, lwgeom2->type=%d", tree1_type, lwgeom2->type); @@ -99,17 +99,9 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2 lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1); lwgeom_free(lwgeom1); } - - /* Need one point from the candidate geometry */ - 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; - } - + /* Flip the candidate point into geographics */ - geographic_point_init(in_point.x, in_point.y, &in_gpoint); + geographic_point_init(in_point->x, in_point->y, &in_gpoint); geog2cart(&in_gpoint, &in_point3d); /* If the candidate isn't in the tree box, it's not in the tree area */ @@ -123,8 +115,8 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2 { POINT2D pt2d_outside; /* latlon */ POINT2D pt2d_inside; - pt2d_inside.x = in_point.x; - pt2d_inside.y = in_point.y; + pt2d_inside.x = in_point->x; + pt2d_inside.y = in_point->y; /* Calculate a definitive outside point */ gbox_pt_outside(&gbox1, &pt2d_outside); POSTGIS_DEBUGF(3, "p2d_inside=POINT(%g %g) p2d_outside=POINT(%g %g)", pt2d_inside.x, pt2d_inside.y, pt2d_outside.x, pt2d_outside.y); @@ -133,81 +125,92 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2 return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL); } } - /* If the un-tree'd argument is a polygon and the tree'd argument isn't, we need to do a */ - /* standard P-i-P on the un-tree'd side. */ - else if ( lwgeom2->type == POLYGONTYPE || lwgeom2->type == MULTIPOLYGONTYPE ) - { - int result; - LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1); - LWGEOM* lwpoint; - POINT4D p4d; - POSTGIS_DEBUG(3, "tree1 not polygonal, but lwgeom2 is, calculating using lwgeom_covers_lwgeom_sphere"); - - if ( LW_FAILURE == lwgeom_startpoint(lwgeom1, &p4d) ) - { - lwgeom_free(lwgeom1); - lwerror("CircTreePIP unable to get lwgeom_startpoint"); - return LW_FALSE; - } - lwpoint = lwpoint_as_lwgeom(lwpoint_make(lwgeom_get_srid(lwgeom1), lwgeom_has_z(lwgeom1), lwgeom_has_m(lwgeom1), &p4d)); - result = lwgeom_covers_lwgeom_sphere(lwgeom2, lwpoint); - lwgeom_free(lwgeom1); - lwgeom_free(lwpoint); - return result; - } else { - POSTGIS_DEBUG(3, "neither tree1 nor lwgeom2 polygonal, so CircTreePIP returning FALSE"); + POSTGIS_DEBUG(3, "tree1 not polygonal, so CircTreePIP returning FALSE"); return LW_FALSE; } } -int -geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance) + +static int +geography_distance_cache_tolerance(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance) { CircTreeGeomCache* tree_cache = NULL; + + int type1 = gserialized_get_type(g1); + int type2 = gserialized_get_type(g2); Assert(distance); /* Two points? Get outa here... */ - if ( (gserialized_get_type(g1) == POINTTYPE) && (gserialized_get_type(g2) == POINTTYPE) ) + if ( type1 == POINTTYPE && type2 == POINTTYPE ) return LW_FAILURE; /* Fetch/build our cache, if appropriate, etc... */ tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2); - if ( tree_cache && tree_cache->argnum && tree_cache->index ) + /* OK, we have an index at the ready! Use it for the one tree argument and */ + /* fill in the other tree argument */ + if ( tree_cache && tree_cache->argnum && tree_cache->index ) { - CIRC_NODE* circ_tree = NULL; - const GSERIALIZED* g = NULL; + CIRC_NODE* circtree_cached = tree_cache->index; + CIRC_NODE* circtree = NULL; + const GSERIALIZED* g_cached; + const GSERIALIZED* g; LWGEOM* lwgeom = NULL; - + int geomtype_cached; + int geomtype; + POINT4D p4d; + /* We need to dynamically build a tree for the uncached side of the function call */ if ( tree_cache->argnum == 1 ) { - lwgeom = lwgeom_from_gserialized(g2); - g = g1; + g_cached = g1; + g = g2; + geomtype_cached = type1; + geomtype = type2; } else if ( tree_cache->argnum == 2 ) { - lwgeom = lwgeom_from_gserialized(g1); - g = g2; + g_cached = g2; + g = g1; + geomtype_cached = type2; + geomtype = type1; } else - lwerror("geography_distance_cache failed! This will never happen!"); - - if ( LW_TRUE == CircTreePIP(tree_cache->index, g, lwgeom) ) { - *distance = 0.0; - lwgeom_free(lwgeom); - return LW_SUCCESS; + lwerror("geography_distance_cache this cannot happen!"); + return LW_FAILURE; } - /* We do tree/tree distance, so turn the candidate geometry into a tree */ - circ_tree = lwgeom_calculate_circ_tree(lwgeom); - /* Calculate tree/tree distance */ - *distance = circ_tree_distance_tree(tree_cache->index, circ_tree, s, FP_TOLERANCE); - circ_tree_free(circ_tree); + lwgeom = lwgeom_from_gserialized(g); + if ( geomtype_cached == POLYGONTYPE || geomtype_cached == MULTIPOLYGONTYPE ) + { + lwgeom_startpoint(lwgeom, &p4d); + if ( CircTreePIP(circtree_cached, g_cached, &p4d) ) + { + *distance = 0.0; + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + } + + circtree = lwgeom_calculate_circ_tree(lwgeom); + if ( geomtype == POLYGONTYPE || geomtype == MULTIPOLYGONTYPE ) + { + circ_tree_get_point(circtree_cached, &p4d); + if ( CircTreePIP(circtree, g, &p4d) ) + { + *distance = 0.0; + circ_tree_free(circtree); + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + } + + *distance = circ_tree_distance_tree(circtree_cached, circtree, s, tolerance); + circ_tree_free(circtree); lwgeom_free(lwgeom); return LW_SUCCESS; } @@ -217,65 +220,27 @@ geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, co } } + +int +geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance) +{ + return geography_distance_cache_tolerance(fcinfo, g1, g2, s, FP_TOLERANCE, distance); +} + int geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin) { - CircTreeGeomCache* tree_cache = NULL; double distance; - - Assert(dwithin); - - /* Two points? Get outa here... */ - if ( (gserialized_get_type(g1) == POINTTYPE) && (gserialized_get_type(g2) == POINTTYPE) ) - return LW_FAILURE; - - /* Fetch/build our cache, if appropriate, etc... */ - tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2); - - if ( tree_cache && tree_cache->argnum && tree_cache->index ) + /* TODO!!! Why does the tolerance stopper in the circ_tree_distance_tree_internal arbitrarily screw up? */ +/* if ( LW_SUCCESS == geography_distance_cache_tolerance(fcinfo, g1, g2, s, tolerance, &distance) ) */ + if ( LW_SUCCESS == geography_distance_cache_tolerance(fcinfo, g1, g2, s, FP_TOLERANCE, &distance) ) { - CIRC_NODE* circ_tree = NULL; - const GSERIALIZED* g = NULL; - LWGEOM* lwgeom = NULL; - - /* We need to dynamically build a tree for the uncached side of the function call */ - if ( tree_cache->argnum == 1 ) - { - lwgeom = lwgeom_from_gserialized(g2); - g = g1; - } - else if ( tree_cache->argnum == 2 ) - { - lwgeom = lwgeom_from_gserialized(g1); - g = g2; - } - else - lwerror("geography_dwithin_cache failed! This will never happen!"); - - if ( LW_TRUE == CircTreePIP(tree_cache->index, g, lwgeom) ) - { - *dwithin = LW_TRUE; - lwgeom_free(lwgeom); - return LW_SUCCESS; - } - - /* We do tree/tree distance, so turn the candidate geometry into a tree */ - circ_tree = lwgeom_calculate_circ_tree(lwgeom); - /* Calculate tree/tree distance */ - distance = circ_tree_distance_tree(tree_cache->index, circ_tree, s, tolerance); *dwithin = (distance <= tolerance ? LW_TRUE : LW_FALSE); - circ_tree_free(circ_tree); - lwgeom_free(lwgeom); return LW_SUCCESS; } - else - { - return LW_FAILURE; - } - + return LW_FAILURE; } - - + int geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance) { @@ -283,13 +248,16 @@ geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHE CIRC_NODE* circ_tree2 = NULL; LWGEOM* lwgeom1 = NULL; LWGEOM* lwgeom2 = NULL; + POINT4D pt1, pt2; lwgeom1 = lwgeom_from_gserialized(g1); lwgeom2 = lwgeom_from_gserialized(g2); circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1); circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2); + lwgeom_startpoint(lwgeom1, &pt1); + lwgeom_startpoint(lwgeom2, &pt2); - if ( CircTreePIP(circ_tree1, g1, lwgeom2) || CircTreePIP(circ_tree2, g2, lwgeom1) ) + if ( CircTreePIP(circ_tree1, g1, &pt2) || CircTreePIP(circ_tree2, g2, &pt1) ) { *distance = 0.0; } diff --git a/regress/tickets.sql b/regress/tickets.sql index a8dccabf1..7ed4977d0 100644 --- a/regress/tickets.sql +++ b/regress/tickets.sql @@ -850,5 +850,14 @@ SELECT '#2168', ST_Distance(g1,g2)::numeric(16,8) As dist_g1_g2, ST_Distance(g FROM (SELECT 'POINT(18.5107234 54.7587757)'::geography As g1, 'POINT(18.58218 54.7344227)'::geography As g2) As a; +-- #2556 -- + +CREATE TABLE images (id integer, name varchar, extent geography(POLYGON,4326)); +INSERT INTO images VALUES (47409, 'TDX-1_2010-10-06T19_44_2375085', 'SRID=4326;POLYGON((-59.4139571913088 82.9486103943668,-57.3528882462655 83.1123152898828,-50.2302874208478 81.3740574826097,-51.977353304689 81.2431047148532,-59.4139571913088 82.9486103943668))'::geography); +INSERT INTO images VALUES (1, 'first_image', 'SRID=4326;POLYGON((-162.211667 88.046667,-151.190278 87.248889,-44.266389 74.887778,-40.793889 75.043333,-162.211667 88.046667))'::geography); +SELECT '#2556' AS ticket, id, round(ST_Distance(extent, 'SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))'::geography)) from images; +DROP TABLE images; + + -- Clean up DELETE FROM spatial_ref_sys; diff --git a/regress/tickets_expected b/regress/tickets_expected index cee55ef1f..56955a799 100644 --- a/regress/tickets_expected +++ b/regress/tickets_expected @@ -253,3 +253,5 @@ ERROR: invalid GML representation #2424|MULTILINESTRING((0 0,10 0,24 3,30 10)) #2427|POINT(-1 0) #2168|5340.76237395|5340.76237395|0 +#2556|47409|20623 +#2556|1|0 -- 2.50.1