From ca02d36dd33cba17638b0ac761e54c9e9f682574 Mon Sep 17 00:00:00 2001 From: Chris Hodgson Date: Thu, 12 May 2011 18:51:57 +0000 Subject: [PATCH] merged fix from r7136 in 1.5 branch, fixes broken point_in_multipolygon_rtree, for #884 git-svn-id: http://svn.osgeo.org/postgis/trunk@7137 b70326c6-7e19-0410-871a-916f4a2858ee --- postgis/lwgeom_functions_analytic.c | 84 ++++++++++++++++------------- postgis/lwgeom_functions_analytic.h | 23 ++++++++ postgis/lwgeom_geos.c | 18 +++---- postgis/lwgeom_rtree.c | 45 ++++++++-------- postgis/lwgeom_rtree.h | 2 +- 5 files changed, 100 insertions(+), 72 deletions(-) create mode 100644 postgis/lwgeom_functions_analytic.h diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c index e0e9e9fdc..17fb44a5c 100644 --- a/postgis/lwgeom_functions_analytic.c +++ b/postgis/lwgeom_functions_analytic.c @@ -17,6 +17,7 @@ #include "math.h" #include "lwgeom_rtree.h" #include "lwalgorithm.h" +#include "lwgeom_functions_analytic.h" /*********************************************************************** @@ -36,11 +37,7 @@ Datum ST_LocateBetweenElevations(PG_FUNCTION_ARGS); double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point); int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point); int point_in_ring(POINTARRAY *pts, POINT2D *point); -int point_in_polygon(LWPOLY *polygon, LWPOINT *point); -int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point); int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point); -int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point); -int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point); PG_FUNCTION_INFO_V1(LWGEOM_simplify2d); @@ -1231,54 +1228,65 @@ int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point) * return 0 if point on boundary * return 1 if point inside polygon * - * Expected **root order is all the exterior rings first, then all the holes - * - * TODO: this could be made slightly more efficient by ordering the rings in - * EIIIEIIIEIEI order (exterior/interior) and including list of exterior ring - * positions on the cache object. + * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI */ -int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point) +int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point) { - int i; + int i, p, r, in_ring; POINT2D pt; int result = -1; - LWDEBUGF(2, "point_in_multipolygon_rtree called for %p %d %d %p.", root, polyCount, ringCount, point); + LWDEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point); getPoint2d_p(point->point, 0, &pt); /* assume bbox short-circuit has already been attempted */ - /* is the point inside (not outside) any of the exterior rings? */ - for ( i = 0; i < polyCount; i++ ) - { - int in_ring = point_in_ring_rtree(root[i], &pt); - LWDEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", i, in_ring); - if ( in_ring != -1 ) /* not outside this ring */ - { - LWDEBUG(3, "point_in_multipolygon_rtree: inside exterior ring."); - result = in_ring; - break; - } - } - - if ( result == -1 ) /* strictly outside all rings */ - return result; + i = 0; /* the current index into the root array */ - /* ok, it's in a ring, but if it's in a hole it's still outside */ - for ( i = polyCount; i < ringCount; i++ ) + /* is the point inside any of the sub-polygons? */ + for ( p = 0; p < polyCount; p++ ) { - int in_ring = point_in_ring_rtree(root[i], &pt); - LWDEBUGF(4, "point_in_multipolygon_rtree: hole (%d), point_in_ring returned %d", i, in_ring); - if ( in_ring == 1 ) /* completely inside hole */ - { - LWDEBUGF(3, "point_in_multipolygon_rtree: within hole %d.", i); - return -1; - } - if ( in_ring == 0 ) /* on the boundary of a hole */ + in_ring = point_in_ring_rtree(root[i], &pt); + LWDEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring); + if ( in_ring == -1 ) /* outside the exterior ring */ { - result = 0; + LWDEBUG(3, "point_in_multipolygon_rtree: outside exterior ring."); + continue; } + if ( in_ring == 0 ) /* on the boundary */ + { + LWDEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p); + return 0; + } + + result = in_ring; + + for(r=1; r outside the polygon */ + { + LWDEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p); + result = -1; + break; + } + if (in_ring == 0) /* on the edge of a hole */ + { + LWDEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p); + return 0; + } + } + /* if we have a positive result, we can short-circuit and return it */ + if ( result != -1) + { + return result; + } + /* increment the index by the total number of rings in the sub-poly */ + /* we do this here in case we short-cutted out of the poly before looking at all the rings */ + i += ringCounts[p]; } + return result; /* -1 = outside, 0 = boundary, 1 = inside */ } diff --git a/postgis/lwgeom_functions_analytic.h b/postgis/lwgeom_functions_analytic.h new file mode 100644 index 000000000..d39bf797b --- /dev/null +++ b/postgis/lwgeom_functions_analytic.h @@ -0,0 +1,23 @@ +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2001-2011 Refractions Research Inc. + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include "lwgeom_rtree.h" + +/* +** Public prototypes for analytic functions. +*/ + +int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point); +int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point); +int point_in_polygon(LWPOLY *polygon, LWPOINT *point); +int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *pont); + diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c index 2e67bd43c..a01571152 100644 --- a/postgis/lwgeom_geos.c +++ b/postgis/lwgeom_geos.c @@ -17,6 +17,7 @@ #include "lwgeom_rtree.h" #include "lwgeom_geos_prepared.h" #include "funcapi.h" +#include "lwgeom_functions_analytic.h" #include #include @@ -72,13 +73,6 @@ Datum pgis_union_geometry_array_old(PG_FUNCTION_ARGS); Datum pgis_union_geometry_array(PG_FUNCTION_ARGS); -/** @todo TODO: move these to a lwgeom_functions_analytic.h - */ -int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point); -int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point); -int point_in_polygon(LWPOLY *polygon, LWPOINT *point); -int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *pont); - /* ** Prototypes end */ @@ -2121,7 +2115,7 @@ Datum contains(PG_FUNCTION_ARGS) if ( poly_cache->ringIndices ) { - result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point); } else if ( type1 == POLYGONTYPE ) { @@ -2364,7 +2358,7 @@ Datum covers(PG_FUNCTION_ARGS) if ( poly_cache->ringIndices ) { - result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point); } else if ( type1 == POLYGONTYPE ) { @@ -2519,7 +2513,7 @@ Datum within(PG_FUNCTION_ARGS) if ( poly_cache->ringIndices ) { - result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point); } else if ( type2 == POLYGONTYPE ) { @@ -2666,7 +2660,7 @@ Datum coveredby(PG_FUNCTION_ARGS) if ( poly_cache->ringIndices ) { - result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point); } else if ( type2 == POLYGONTYPE ) { @@ -2901,7 +2895,7 @@ Datum intersects(PG_FUNCTION_ARGS) if ( poly_cache->ringIndices ) { - result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point); } else if ( polytype == POLYGONTYPE ) { diff --git a/postgis/lwgeom_rtree.c b/postgis/lwgeom_rtree.c index a9da98644..63fc3443f 100644 --- a/postgis/lwgeom_rtree.c +++ b/postgis/lwgeom_rtree.c @@ -212,17 +212,23 @@ void freeTree(RTREE_NODE *root) */ void clearCache(RTREE_POLY_CACHE *cache) { - int i; + int g, r, i; LWDEBUGF(2, "clearCache called for %p", cache); - for (i = 0; i < cache->ringCount; i++) - { - freeTree(cache->ringIndices[i]); + i = 0; + for (g = 0; g < cache->polyCount; g++) + { + for (r = 0; r < cache->ringCounts; r++) + { + freeTree(cache->ringIndices[i]); + i++; + } } lwfree(cache->ringIndices); + lwfree(cache->ringCounts); lwfree(cache->poly); cache->poly = 0; cache->ringIndices = 0; - cache->ringCount = 0; + cache->ringCounts = 0; cache->polyCount = 0; } @@ -391,7 +397,7 @@ RTREE_POLY_CACHE * createCache() RTREE_POLY_CACHE *result; result = lwalloc(sizeof(RTREE_POLY_CACHE)); result->polyCount = 0; - result->ringCount = 0; + result->ringCounts = 0; result->ringIndices = 0; result->poly = 0; result->type = 1; @@ -400,7 +406,7 @@ RTREE_POLY_CACHE * createCache() void populateCache(RTREE_POLY_CACHE *currentCache, LWGEOM *lwgeom, PG_LWGEOM *serializedPoly) { - int i, j, k, length; + int i, p, r, length; LWMPOLY *mpoly; LWPOLY *poly; int nrings; @@ -415,28 +421,24 @@ void populateCache(RTREE_POLY_CACHE *currentCache, LWGEOM *lwgeom, PG_LWGEOM *se /* ** Count the total number of rings. */ + currentCache->polyCount = mpoly->ngeoms; + currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms); for ( i = 0; i < mpoly->ngeoms; i++ ) { + currentCache->ringCounts[i] = mpoly->geoms[i]->nrings; nrings += mpoly->geoms[i]->nrings; } - currentCache->polyCount = mpoly->ngeoms; - currentCache->ringCount = nrings; currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings); /* - ** Load the exterior rings onto the ringIndices array first - */ - for ( i = 0; i < mpoly->ngeoms; i++ ) - { - currentCache->ringIndices[i] = createTree(mpoly->geoms[i]->rings[0]); - } - /* - ** Load the interior rings (holes) onto ringIndices next + ** Load the array in geometry order, each outer ring followed by the inner rings + ** associated with that outer ring */ - for ( j = 0; j < mpoly->ngeoms; j++ ) + i = 0; + for ( p = 0; p < mpoly->ngeoms; p++ ) { - for ( k = 1; k < mpoly->geoms[j]->nrings; k++ ) + for ( r = 0; r < mpoly->geoms[p]->nrings; r++ ) { - currentCache->ringIndices[i] = createTree(mpoly->geoms[j]->rings[k]); + currentCache->ringIndices[i] = createTree(mpoly->geoms[p]->rings[r]); i++; } } @@ -446,7 +448,8 @@ void populateCache(RTREE_POLY_CACHE *currentCache, LWGEOM *lwgeom, PG_LWGEOM *se LWDEBUG(2, "populateCache POLYGON"); poly = (LWPOLY *)lwgeom; currentCache->polyCount = 1; - currentCache->ringCount = poly->nrings; + currentCache->ringCounts = lwalloc(sizeof(int)); + currentCache->ringCounts[0] = poly->nrings; /* ** Just load the rings on in order */ diff --git a/postgis/lwgeom_rtree.h b/postgis/lwgeom_rtree.h index c6469cbe8..6933ac383 100644 --- a/postgis/lwgeom_rtree.h +++ b/postgis/lwgeom_rtree.h @@ -49,7 +49,7 @@ typedef struct { char type; RTREE_NODE **ringIndices; - int ringCount; + int* ringCounts; int polyCount; PG_LWGEOM *poly; } -- 2.50.1