From 1e2efa6144d9d9a9c1db28d98a8c8c1618f58c01 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Mon, 29 Sep 2008 18:25:00 +0000 Subject: [PATCH] Backport memory improvements and multi-polygon support for p-i-p shortcuts from trunk. git-svn-id: http://svn.osgeo.org/postgis/branches/1.3@3028 b70326c6-7e19-0410-871a-916f4a2858ee --- lwgeom/lwgeom_functions_analytic.c | 232 +++++++++------ lwgeom/lwgeom_functions_lrs.c | 15 +- lwgeom/lwgeom_geos_c.c | 456 +++++++++++++---------------- lwgeom/lwgeom_rtree.c | 196 +++++++++---- lwgeom/lwgeom_rtree.h | 7 +- 5 files changed, 489 insertions(+), 417 deletions(-) diff --git a/lwgeom/lwgeom_functions_analytic.c b/lwgeom/lwgeom_functions_analytic.c index a88e675a0..a1dfa128c 100644 --- a/lwgeom/lwgeom_functions_analytic.c +++ b/lwgeom/lwgeom_functions_analytic.c @@ -1136,7 +1136,7 @@ int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point) * return 1 iff point is inside ring pts * return 0 iff point is on ring pts */ -int point_in_ring(RTREE_NODE *root, POINT2D *point) +int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point) { int wn = 0; int i; @@ -1224,7 +1224,7 @@ int point_in_ring(RTREE_NODE *root, POINT2D *point) * return 1 iff point is inside ring pts * return 0 iff point is on ring pts */ -int point_in_ring_deprecated(POINTARRAY *pts, POINT2D *point) +int point_in_ring(POINTARRAY *pts, POINT2D *point) { int wn = 0; int i; @@ -1306,7 +1306,7 @@ int point_in_ring_deprecated(POINTARRAY *pts, POINT2D *point) * return 0 iff point outside polygon or on boundary * return 1 iff point inside polygon */ -int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) +int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point) { int i; POINT2D pt; @@ -1318,7 +1318,7 @@ int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) getPoint2d_p(point->point, 0, &pt); /* assume bbox short-circuit has already been attempted */ - if(point_in_ring(root[0], &pt) != 1) + if(point_in_ring_rtree(root[0], &pt) != 1) { #ifdef PGIS_DEBUG lwnotice("point_in_polygon: outside exterior ring."); @@ -1328,7 +1328,7 @@ int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) for(i=1; ipoint, 0, &pt); - /* assume bbox short-circuit has already been attempted */ - - ring = polygon->rings[0]; - /* root = createTree(ring); */ - /* if(point_in_ring(root, &pt) != 1) */ - if(point_in_ring_deprecated(polygon->rings[0], &pt) != 1) - { -#ifdef PGIS_DEBUG - lwnotice("point_in_polygon: outside exterior ring."); -#endif - return 0; - } + getPoint2d_p(point->point, 0, &pt); + /* assume bbox short-circuit has already been attempted */ - for(i=1; inrings; i++) - { - ring = polygon->rings[i]; - /* root = createTree(ring); */ - /* if(point_in_ring(root, &pt) != -1) */ - if(point_in_ring_deprecated(polygon->rings[i], &pt) != -1) - { -#ifdef PGIS_DEBUG - lwnotice("point_in_polygon: within hole %d.", i); + /* 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); +#ifdef PGIS_DEBUG_CALLS + lwnotice("point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", i, in_ring); #endif - return 0; - } - } - return 1; + if( in_ring != -1 ) /* not outside this ring */ + { +#ifdef PGIS_DEBUG_CALLS + lwnotice("point_in_multipolygon_rtree: inside exterior ring."); +#endif + result = in_ring; + break; + } + } + + if( result == -1 ) /* strictly outside all rings */ + return result; + + /* ok, it's in a ring, but if it's in a hole it's still outside */ + for( i = polyCount; i < ringCount; i++ ) + { + int in_ring = point_in_ring_rtree(root[i], &pt); +#ifdef PGIS_DEBUG_CALLS + lwnotice("point_in_multipolygon_rtree: hole (%d), point_in_ring returned %d", i, in_ring); +#endif + if( in_ring == 1 ) /* completely inside hole */ + { +#ifdef PGIS_DEBUG_CALLS + lwnotice("point_in_multipolygon_rtree: within hole %d.", i); +#endif + return -1; + } + if( in_ring == 0 ) /* on the boundary of a hole */ + { + result = 0; + } + } + return result; /* -1 = outside, 0 = boundary, 1 = inside */ + } + + /* - * return 0 iff point inside polygon or on boundary - * return 1 iff point outside polygon + * return -1 iff point outside polygon + * return 0 iff point on boundary + * return 1 iff point inside polygon */ -int point_outside_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) +int point_in_polygon(LWPOLY *polygon, LWPOINT *point) { - int i; + int i, result, in_ring; + POINTARRAY *ring; POINT2D pt; #ifdef PGIS_DEBUG_CALLS - lwnotice("point_outside_polygon called."); + lwnotice("point_in_polygon_deprecated called."); #endif getPoint2d_p(point->point, 0, &pt); /* assume bbox short-circuit has already been attempted */ - if(point_in_ring(root[0], &pt) == -1) + ring = polygon->rings[0]; + in_ring = point_in_ring(polygon->rings[0], &pt); + if( in_ring == -1) /* outside the exterior ring */ { #ifdef PGIS_DEBUG - lwnotice("point_outside_polygon: outside exterior ring."); + lwnotice("point_in_polygon: outside exterior ring."); #endif - return 1; + return -1; } + result = in_ring; - for(i=1; inrings; i++) { - if(point_in_ring(root[i], &pt) == 1) + ring = polygon->rings[i]; + in_ring = point_in_ring(polygon->rings[i], &pt); + if(in_ring == 1) /* inside a hole => outside the polygon */ { #ifdef PGIS_DEBUG - lwnotice("point_outside_polygon: within hole %d.", i); + lwnotice("point_in_polygon: within hole %d.", i); #endif - return 1; + return -1; + } + if(in_ring == 0) /* on the edge of a hole */ + { +#ifdef PGIS_DEBUG + lwnotice("point_in_polygon: on edge of hole %d.", i); +#endif + return 0; } } - return 0; + return result; /* -1 = outside, 0 = boundary, 1 = inside */ } /* - * return 0 iff point inside polygon or on boundary - * return 1 iff point outside polygon + * return -1 iff point outside multipolygon + * return 0 iff point on multipolygon boundary + * return 1 iff point inside multipolygon */ -int point_outside_polygon_deprecated(LWPOLY *polygon, LWPOINT *point) +int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point) { - int i; + int i, j, result, in_ring; POINTARRAY *ring; POINT2D pt; -#ifdef PGIS_DEBUG_CALLS - lwnotice("point_outside_polygon_deprecated called."); +#ifdef PGIS_DEBUG + lwnotice("point_in_polygon called."); #endif getPoint2d_p(point->point, 0, &pt); /* assume bbox short-circuit has already been attempted */ - - ring = polygon->rings[0]; - /* root = createTree(ring); */ - /* if(point_in_ring(root, &pt) == -1) */ - if(point_in_ring_deprecated(ring, &pt) == -1) - { + + result = -1; + + for(j = 0; j < mpolygon->ngeoms; j++ ) + { + + LWPOLY *polygon = mpolygon->geoms[j]; + ring = polygon->rings[0]; + in_ring = point_in_ring(polygon->rings[0], &pt); + if( in_ring == -1) /* outside the exterior ring */ + { #ifdef PGIS_DEBUG - lwnotice("point_outside_polygon_deprecated: outside exterior ring."); + lwnotice("point_in_polygon: outside exterior ring."); #endif - return 1; - } + continue; + } + if( in_ring == 0 ) + { + return 0; + } - for(i=1; inrings; i++) - { + result = in_ring; + + for(i=1; inrings; i++) + { ring = polygon->rings[i]; - /* root = createTree(ring); */ - /* if(point_in_ring(root, &pt) == 1) */ - if(point_in_ring_deprecated(ring, &pt) == 1) + in_ring = point_in_ring(polygon->rings[i], &pt); + if(in_ring == 1) /* inside a hole => outside the polygon */ { #ifdef PGIS_DEBUG - lwnotice("point_outside_polygon_deprecated: within hole %d.", i); + lwnotice("point_in_polygon: within hole %d.", i); #endif - return 1; + result = -1; + break; } - } - return 0; -} - - -/* - * return 0 iff point is outside every polygon - */ -/* Not yet functional. -int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point) -{ - int i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("point_in_multipolygon called."); + if(in_ring == 0) /* on the edge of a hole */ + { +#ifdef PGIS_DEBUG + lwnotice("point_in_polygon: on edge of hole %d.", i); #endif - - for(i=1; ingeoms; i++) - { - if(point_in_polygon((LWPOLY *)mpolygon->geoms[i], point)!=0) return 1; - } - return 0; + return 0; + } + } + if( result != -1) + { + return result; + } + } + return result; } -*/ /******************************************************************************* diff --git a/lwgeom/lwgeom_functions_lrs.c b/lwgeom/lwgeom_functions_lrs.c index f7e1e8a53..dc9f62ce8 100644 --- a/lwgeom/lwgeom_functions_lrs.c +++ b/lwgeom/lwgeom_functions_lrs.c @@ -16,8 +16,8 @@ #include "lwgeom_pg.h" #include "math.h" -#define DEBUG_LRS 0 -#define DEBUG_INTERPOLATION 0 +#define DEBUG_LRS 1 +#define DEBUG_INTERPOLATION 1 Datum LWGEOM_locate_between_m(PG_FUNCTION_ARGS); @@ -151,7 +151,7 @@ clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1) * numbers. * * The two points must be equal anyway. - */ + if ( m0 == m1 ) { memcpy(p2, p1, sizeof(POINT4D)); @@ -159,6 +159,7 @@ clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1) else ret |= 0x0100; return ret; } +XXX */ /* * Second point out of range, project @@ -338,7 +339,7 @@ lwpoint_locate_between_m(LWPOINT *lwpoint, double m0, double m1) */ static LWGEOM * lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1) -{ +{ POINTARRAY *ipa=lwline_in->points; int i; LWGEOM **geoms; @@ -377,6 +378,7 @@ lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1) /* This is a point */ if ( pa->npoints == 1 ) { + lwpoint=lwalloc(sizeof(LWPOINT)); lwpoint->type=lwgeom_makeType_full( TYPE_HASZ(pa->dims), @@ -414,6 +416,8 @@ lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1) lwerror("ptarray_locate_between_m returned a POINARRAY set containing POINTARRAY with 0 points"); } + + } if ( ngeoms == 1 ) @@ -427,6 +431,9 @@ lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1) else if ( typeflag == 2 ) outtype=MULTILINETYPE; else outtype = COLLECTIONTYPE; + lwnotice(" XXX lwline_locate_between_m: %s", lwgeom_to_ewkt((LWGEOM *)lwcollection_construct(outtype, + lwline_in->SRID, NULL, ngeoms, geoms))); + return (LWGEOM *)lwcollection_construct(outtype, lwline_in->SRID, NULL, ngeoms, geoms); } diff --git a/lwgeom/lwgeom_geos_c.c b/lwgeom/lwgeom_geos_c.c index 440b4e009..f74247e91 100644 --- a/lwgeom/lwgeom_geos_c.c +++ b/lwgeom/lwgeom_geos_c.c @@ -46,6 +46,8 @@ /* #define UNITE_USING_BUFFER 1 */ /* #define MAXGEOMSPOINTS 21760 */ +/* PROTOTYPES start */ + Datum relate_full(PG_FUNCTION_ARGS); Datum relate_pattern(PG_FUNCTION_ARGS); Datum disjoint(PG_FUNCTION_ARGS); @@ -85,6 +87,12 @@ GEOSGeom POSTGIS2GEOS(PG_LWGEOM *g); GEOSGeom LWGEOM2GEOS(LWGEOM *g); void errorIfGeometryCollection(PG_LWGEOM *g1, PG_LWGEOM *g2); +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 */ PG_FUNCTION_INFO_V1(postgis_geos_version); Datum postgis_geos_version(PG_FUNCTION_ARGS) @@ -1430,12 +1438,6 @@ Datum overlaps(PG_FUNCTION_ARGS) PG_RETURN_BOOL(result); } -int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point); -int point_in_polygon_deprecated(LWPOLY *polygon, LWPOINT *point); -int point_outside_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point); -int point_outside_polygon_deprecated(LWPOLY *polygon, LWPOINT *point); -int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point); - PG_FUNCTION_INFO_V1(contains); Datum contains(PG_FUNCTION_ARGS) @@ -1446,7 +1448,7 @@ Datum contains(PG_FUNCTION_ARGS) bool result; BOX2DFLOAT4 box1, box2; int type1, type2; - LWPOLY *poly; + LWGEOM *lwgeom; /* LWMPOLY *mpoly; */ LWPOINT *point; RTREE_POLY_CACHE *poly_cache; @@ -1467,13 +1469,14 @@ Datum contains(PG_FUNCTION_ARGS) * geom1 bounding box we can prematurely return FALSE. * Do the test IFF BOUNDING BOX AVAILABLE. */ - if ( getbox2d_p(SERIALIZED_FORM(geom1), &box1) && - getbox2d_p(SERIALIZED_FORM(geom2), &box2) ) + if ( getbox2d_p(SERIALIZED_FORM(geom1), &box1) && + getbox2d_p(SERIALIZED_FORM(geom2), &box2) ) { - if ( box2.xmin < box1.xmin ) PG_RETURN_BOOL(FALSE); - if ( box2.xmax > box1.xmax ) PG_RETURN_BOOL(FALSE); - if ( box2.ymin < box1.ymin ) PG_RETURN_BOOL(FALSE); - if ( box2.ymax > box1.ymax ) PG_RETURN_BOOL(FALSE); + if ( ( box2.xmin < box1.xmin ) || ( box2.xmax > box1.xmax ) || + ( box2.ymin < box1.ymin ) || ( box2.ymax > box1.ymax ) ) + { + PG_RETURN_BOOL(FALSE); + } } /* * short-circuit 2: if geom2 is a point and geom1 is a polygon @@ -1481,13 +1484,13 @@ Datum contains(PG_FUNCTION_ARGS) */ type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POLYGONTYPE && type2 == POINTTYPE) - { + if((type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE) && type2 == POINTTYPE) + { #ifdef PGIS_DEBUG lwnotice("Point in Polygon test requested...short-circuiting."); #endif - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); #ifdef PGIS_DEBUG lwnotice("Precall point_in_polygon %p, %p", poly, point); @@ -1499,53 +1502,40 @@ Datum contains(PG_FUNCTION_ARGS) * future use, then switch back to the local context. */ old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom1), fcinfo->flinfo->fn_extra); + poly_cache = retrieveCache(lwgeom, SERIALIZED_FORM(geom1), fcinfo->flinfo->fn_extra); fcinfo->flinfo->fn_extra = poly_cache; MemoryContextSwitchTo(old_context); - if(point_in_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } + if( poly_cache->ringIndices ) + { + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + } + else if ( type1 == POLYGONTYPE ) + { + result = point_in_polygon((LWPOLY*)lwgeom, point); + } + else if ( type1 == MULTIPOLYGONTYPE ) + { + result = point_in_multipolygon((LWMPOLY*)lwgeom, point); + } + else { + /* Gulp! Should not be here... */ + elog(ERROR,"Type isn't poly or multipoly!"); + PG_RETURN_NULL(); + } + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + lwgeom_release((LWGEOM *)lwgeom); + lwgeom_release((LWGEOM *)point); + if( result == 1 ) /* completely inside */ + { + PG_RETURN_BOOL(TRUE); + } + else + { + PG_RETURN_BOOL(FALSE); + } } - /* Not yet functional - else if(type1 == MULTIPOLYGONTYPE && type2 == POINTTYPE) - { -#ifdef PGIS_DEBUG - lwnotice("Point in MultiPolygon test requested...short-circuiting."); -#endif - mpoly = lwmpoly_deserialize(SERIALIZED_FORM(geom1)); - point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); - if(point_in_multipolygon(mpoly, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)mpoly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)mpoly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } - } - */ else { #ifdef PGIS_DEBUG @@ -1612,7 +1602,7 @@ Datum covers(PG_FUNCTION_ARGS) bool result; BOX2DFLOAT4 box1, box2; int type1, type2; - LWPOLY *poly; + LWGEOM *lwgeom; /* LWMPOLY *mpoly; */ LWPOINT *point; RTREE_POLY_CACHE *poly_cache; @@ -1648,13 +1638,13 @@ Datum covers(PG_FUNCTION_ARGS) */ type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POLYGONTYPE && type2 == POINTTYPE) + if((type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE) && type2 == POINTTYPE) { #ifdef PGIS_DEBUG lwnotice("Point in Polygon test requested...short-circuiting."); #endif - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); #ifdef PGIS_DEBUG lwnotice("Precall point_in_polygon %p, %p", poly, point); @@ -1666,26 +1656,40 @@ Datum covers(PG_FUNCTION_ARGS) * future use, then switch back to the local context. */ old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom1), fcinfo->flinfo->fn_extra); + poly_cache = retrieveCache(lwgeom, SERIALIZED_FORM(geom1), fcinfo->flinfo->fn_extra); fcinfo->flinfo->fn_extra = poly_cache; MemoryContextSwitchTo(old_context); - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } + if( poly_cache->ringIndices ) + { + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + } + else if ( type1 == POLYGONTYPE ) + { + result = point_in_polygon((LWPOLY*)lwgeom, point); + } + else if ( type1 == MULTIPOLYGONTYPE ) + { + result = point_in_multipolygon((LWMPOLY*)lwgeom, point); + } + else { + /* Gulp! Should not be here... */ + elog(ERROR,"Type isn't poly or multipoly!"); + PG_RETURN_NULL(); + } + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + lwgeom_release((LWGEOM *)lwgeom); + lwgeom_release((LWGEOM *)point); + if( result != -1 ) /* not outside */ + { + PG_RETURN_BOOL(TRUE); + } + else + { + PG_RETURN_BOOL(FALSE); + } } else { @@ -1747,7 +1751,7 @@ Datum within(PG_FUNCTION_ARGS) GEOSGeom g1,g2; bool result; BOX2DFLOAT4 box1, box2; - LWPOLY *poly; + LWGEOM *lwgeom; LWPOINT *point; int type1, type2; MemoryContext old_context; @@ -1782,14 +1786,14 @@ Datum within(PG_FUNCTION_ARGS) */ type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POINTTYPE && type2 == POLYGONTYPE) + if((type2 == POLYGONTYPE || type2 == MULTIPOLYGONTYPE) && type1 == POINTTYPE) { #ifdef PGIS_DEBUG lwnotice("Point in Polygon test requested...short-circuiting."); #endif point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom2)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom2)); /* * Switch the context to the function-scope context, @@ -1797,26 +1801,40 @@ Datum within(PG_FUNCTION_ARGS) * future use, then switch back to the local context. */ old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); + poly_cache = retrieveCache(lwgeom, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); fcinfo->flinfo->fn_extra = poly_cache; MemoryContextSwitchTo(old_context); - if(point_in_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } + if( poly_cache->ringIndices ) + { + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + } + else if ( type2 == POLYGONTYPE ) + { + result = point_in_polygon((LWPOLY*)lwgeom, point); + } + else if ( type2 == MULTIPOLYGONTYPE ) + { + result = point_in_multipolygon((LWMPOLY*)lwgeom, point); + } + else { + /* Gulp! Should not be here... */ + elog(ERROR,"Type isn't poly or multipoly!"); + PG_RETURN_NULL(); + } + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + lwgeom_release((LWGEOM *)lwgeom); + lwgeom_release((LWGEOM *)point); + if( result == 1 ) /* completely inside */ + { + PG_RETURN_BOOL(TRUE); + } + else + { + PG_RETURN_BOOL(FALSE); + } } initGEOS(lwnotice, lwnotice); @@ -1877,7 +1895,7 @@ Datum coveredby(PG_FUNCTION_ARGS) GEOSGeom g1,g2; bool result; BOX2DFLOAT4 box1, box2; - LWPOLY *poly; + LWGEOM *lwgeom; LWPOINT *point; int type1, type2; MemoryContext old_context; @@ -1917,14 +1935,14 @@ Datum coveredby(PG_FUNCTION_ARGS) */ type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POINTTYPE && type2 == POLYGONTYPE) + if((type2 == POLYGONTYPE || type2 == MULTIPOLYGONTYPE) && type1 == POINTTYPE) { #ifdef PGIS_DEBUG lwnotice("Point in Polygon test requested...short-circuiting."); #endif point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom2)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom2)); /* * Switch the context to the function-scope context, @@ -1932,26 +1950,40 @@ Datum coveredby(PG_FUNCTION_ARGS) * future use, then switch back to the local context. */ old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); + poly_cache = retrieveCache(lwgeom, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); fcinfo->flinfo->fn_extra = poly_cache; MemoryContextSwitchTo(old_context); - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } + if( poly_cache->ringIndices ) + { + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + } + else if ( type2 == POLYGONTYPE ) + { + result = point_in_polygon((LWPOLY*)lwgeom, point); + } + else if ( type2 == MULTIPOLYGONTYPE ) + { + result = point_in_multipolygon((LWMPOLY*)lwgeom, point); + } + else { + /* Gulp! Should not be here... */ + elog(ERROR,"Type isn't poly or multipoly!"); + PG_RETURN_NULL(); + } + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + lwgeom_release((LWGEOM *)lwgeom); + lwgeom_release((LWGEOM *)point); + if( result != -1 ) /* not outside */ + { + PG_RETURN_BOOL(TRUE); + } + else + { + PG_RETURN_BOOL(FALSE); + } } initGEOS(lwnotice, lwnotice); @@ -2085,13 +2117,14 @@ PG_FUNCTION_INFO_V1(intersects); Datum intersects(PG_FUNCTION_ARGS) { PG_LWGEOM *geom1; - PG_LWGEOM *geom2; + PG_LWGEOM *geom2; + uchar *serialized_poly; GEOSGeom g1,g2; bool result; BOX2DFLOAT4 box1, box2; - int type1, type2; + int type1, type2, polytype; LWPOINT *point; - LWPOLY *poly; + LWGEOM *lwgeom; MemoryContext old_context; RTREE_POLY_CACHE *poly_cache; @@ -2125,76 +2158,63 @@ Datum intersects(PG_FUNCTION_ARGS) */ type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POINTTYPE && type2 == POLYGONTYPE) + if( (type1 == POINTTYPE && (type2 == POLYGONTYPE || type2 == MULTIPOLYGONTYPE)) || + (type2 == POINTTYPE && (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE))) { #ifdef PGIS_DEBUG lwnotice("Point in Polygon test requested...short-circuiting."); #endif - point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom2)); - - /* - * Switch the context to the function-scope context, - * retrieve the appropriate cache object, cache it for - * future use, then switch back to the local context. - */ - old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); - fcinfo->flinfo->fn_extra = poly_cache; - MemoryContextSwitchTo(old_context); + if( type1 == POINTTYPE ) { + point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom2)); + serialized_poly = SERIALIZED_FORM(geom2); + polytype = type2; + } else { + point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); + serialized_poly = SERIALIZED_FORM(geom1); + polytype = type1; + } - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) + /* + * Switch the context to the function-scope context, + * retrieve the appropriate cache object, cache it for + * future use, then switch back to the local context. + */ + old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); + poly_cache = retrieveCache(lwgeom, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); + fcinfo->flinfo->fn_extra = poly_cache; + MemoryContextSwitchTo(old_context); + + if( poly_cache->ringIndices ) { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); + result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); } - else + else if ( polytype == POLYGONTYPE ) { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); + result = point_in_polygon((LWPOLY*)lwgeom, point); + } + else if ( polytype == MULTIPOLYGONTYPE ) + { + result = point_in_multipolygon((LWMPOLY*)lwgeom, point); + } + else { + /* Gulp! Should not be here... */ + elog(ERROR,"Type isn't poly or multipoly!"); + PG_RETURN_NULL(); } - } - else if(type1 == POLYGONTYPE && type2 == POINTTYPE) - { - #ifdef PGIS_DEBUG - lwnotice("Point in Polygon test requested...short-circuiting."); -#endif - point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); - - /* - * Switch the context to the function-scope context, - * retrieve the appropriate cache object, cache it for - * future use, then switch back to the local context. - */ - old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom1), fcinfo->flinfo->fn_extra); - fcinfo->flinfo->fn_extra = poly_cache; - MemoryContextSwitchTo(old_context); - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + lwgeom_release((LWGEOM *)lwgeom); + lwgeom_release((LWGEOM *)point); + if( result != -1 ) /* not outside */ { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); PG_RETURN_BOOL(TRUE); } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 1); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); + else { PG_RETURN_BOOL(FALSE); } - } @@ -2329,11 +2349,6 @@ Datum disjoint(PG_FUNCTION_ARGS) GEOSGeom g1,g2; bool result; BOX2DFLOAT4 box1, box2; - int type1, type2; - LWPOLY *poly; - LWPOINT *point; - MemoryContext old_context; - RTREE_POLY_CACHE *poly_cache; #ifdef PROFILE profstart(PROF_QRUN); @@ -2359,83 +2374,6 @@ Datum disjoint(PG_FUNCTION_ARGS) if ( box2.ymin > box1.ymax ) PG_RETURN_BOOL(TRUE); } - /* - * short-circuit 2: if the geoms are a point and a polygon, - * call the point_outside_polygon function. - */ - type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); - type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POINTTYPE && type2 == POLYGONTYPE) - { -#ifdef PGIS_DEBUG - lwnotice("Point outside Polygon test requested...short-circuiting."); -#endif - point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom2)); - - /* - * Switch the context to the function-scope context, - * retrieve the appropriate cache object, cache it for - * future use, then switch back to the local context. - */ - old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); - fcinfo->flinfo->fn_extra = poly_cache; - MemoryContextSwitchTo(old_context); - - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } - } - else if(type1 == POLYGONTYPE && type2 == POINTTYPE) - { -#ifdef PGIS_DEBUG - lwnotice("Point outside Polygon test requested...short-circuiting."); -#endif - point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); - - /* - * Switch the context to the function-scope context, - * retrieve the appropriate cache object, cache it for - * future use, then switch back to the local context. - */ - old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); - poly_cache = retrieveCache(poly, SERIALIZED_FORM(geom2), fcinfo->flinfo->fn_extra); - fcinfo->flinfo->fn_extra = poly_cache; - MemoryContextSwitchTo(old_context); - - if(point_outside_polygon(poly_cache->ringIndices, poly_cache->ringCount, point) == 0) - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(FALSE); - } - else - { - PG_FREE_IF_COPY(geom1, 0); - PG_FREE_IF_COPY(geom2, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)point); - PG_RETURN_BOOL(TRUE); - } - } - initGEOS(lwnotice, lwnotice); #ifdef PROFILE diff --git a/lwgeom/lwgeom_rtree.c b/lwgeom/lwgeom_rtree.c index 66d294b9f..74b8e5325 100644 --- a/lwgeom/lwgeom_rtree.c +++ b/lwgeom/lwgeom_rtree.c @@ -77,6 +77,7 @@ RTREE_NODE *createTree(POINTARRAY *pointArray) } root = nodes[0]; + lwfree(nodes); #ifdef PGIS_DEBUG lwnotice("createTree returning %p", root); @@ -207,16 +208,41 @@ void freeTree(RTREE_NODE *root) #ifdef PGIS_DEBUG_CALLS lwnotice("freeTree called for %p", root); #endif - if(root->leftNode) - freeTree(root->leftNode); - if(root->rightNode) - freeTree(root->rightNode); - lwfree(root->interval); - if(root->segment) - lwfree(root->segment); - lwfree(root); + if(root->leftNode) + freeTree(root->leftNode); + if(root->rightNode) + freeTree(root->rightNode); + lwfree(root->interval); + if(root->segment) { + lwfree(root->segment->points->serialized_pointlist); + lwfree(root->segment->points); + lwgeom_release((LWGEOM *)root->segment); + } + lwfree(root); } +/* + * Free the cache object and all the sub-objects properly. + */ +void clearCache(RTREE_POLY_CACHE *cache) +{ + int i; +#ifdef PGIS_DEBUG_CALLS + lwnotice("clearCache called for %p", cache); +#endif + for(i = 0; i < cache->ringCount; i++) + { + freeTree(cache->ringIndices[i]); + } + lwfree(cache->ringIndices); + lwfree(cache->poly); + cache->poly = 0; + cache->ringIndices = 0; + cache->ringCount = 0; + cache->polyCount = 0; +} + + /* * Retrieves a collection of line segments given the root and crossing value. * The collection is a multilinestring consisting of two point lines @@ -387,28 +413,98 @@ Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS) } -RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly) +RTREE_POLY_CACHE * createCache() { - RTREE_POLY_CACHE *result; - int i, length; + RTREE_POLY_CACHE *result; + result = lwalloc(sizeof(RTREE_POLY_CACHE)); + result->polyCount = 0; + result->ringCount = 0; + result->ringIndices = 0; + result->poly = 0; + return result; +} +void populateCache(RTREE_POLY_CACHE *currentCache, LWGEOM *lwgeom, uchar *serializedPoly) +{ + int i, j, k, length; + LWMPOLY *mpoly; + LWPOLY *poly; + int nrings; + #ifdef PGIS_DEBUG_CALLS - lwnotice("createNewCache called with %p", poly); + lwnotice("populateCache called with cache %p geom %p", currentCache, lwgeom); #endif - result = lwalloc(sizeof(RTREE_POLY_CACHE)); - result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings); - result->ringCount = poly->nrings; - length = lwgeom_size_poly(serializedPoly); - result->poly = lwalloc(length); - memcpy(result->poly, serializedPoly, length); - for(i = 0; i < result->ringCount; i++) - { - result->ringIndices[i] = createTree(poly->rings[i]); - } -#ifdef PGIS_DEBUG - lwnotice("createNewCache returning %p", result); + + if(TYPE_GETTYPE(lwgeom->type) == MULTIPOLYGONTYPE) + { +#ifdef PGIS_DEBUG_CALLS + lwnotice("populateCache MULTIPOLYGON"); +#endif + mpoly = (LWMPOLY *)lwgeom; + nrings = 0; + /* + ** Count the total number of rings. + */ + for( i = 0; i < mpoly->ngeoms; i++ ) + { + 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 + */ + for( j = 0; j < mpoly->ngeoms; j++ ) + { + for( k = 1; k < mpoly->geoms[j]->nrings; k++ ) + { + currentCache->ringIndices[i] = createTree(mpoly->geoms[j]->rings[k]); + i++; + } + } + } + else if ( TYPE_GETTYPE(lwgeom->type) == POLYGONTYPE ) + { +#ifdef PGIS_DEBUG_CALLS + lwnotice("populateCache POLYGON"); +#endif + poly = (LWPOLY *)lwgeom; + currentCache->polyCount = 1; + currentCache->ringCount = poly->nrings; + /* + ** Just load the rings on in order + */ + currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings); + for( i = 0; i < poly->nrings; i++ ) + { + currentCache->ringIndices[i] = createTree(poly->rings[i]); + } + } + else + { + /* Uh oh, shouldn't be here. */ + return; + } + + /* + ** Copy the serialized form of the polygon into the cache so + ** we can test for equality against subsequent polygons. + */ + length = lwgeom_size(serializedPoly); + currentCache->poly = lwalloc(length); + memcpy(currentCache->poly, serializedPoly, length); +#ifdef PGIS_DEBUG_CALLS + lwnotice("populateCache returning %p", currentCache); #endif - return result; + } /* @@ -418,10 +514,10 @@ RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly) * method. The method will allocate memory for the cache it creates, * as well as freeing the memory of any cache that is no longer applicable. */ -RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, +RTREE_POLY_CACHE *retrieveCache(LWGEOM *lwgeom, uchar *serializedPoly, RTREE_POLY_CACHE *currentCache) { - int i, length; + int length; #ifdef PGIS_DEBUG_CALLS lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache); @@ -431,51 +527,35 @@ RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, #ifdef PGIS_DEBUG lwnotice("No existing cache, create one."); #endif - return createNewCache(poly, serializedPoly); + return createCache(); } if(!(currentCache->poly)) { #ifdef PGIS_DEBUG - lwnotice("Cache contains no polygon, creating new cache."); + lwnotice("Cache contains no polygon, populating it."); #endif - return createNewCache(poly, serializedPoly); + populateCache(currentCache, lwgeom, serializedPoly); + return currentCache; } length = lwgeom_size_poly(serializedPoly); - if(lwgeom_size_poly(currentCache->poly) != length) + if(lwgeom_size(currentCache->poly) != length) { #ifdef PGIS_DEBUG lwnotice("Polygon size mismatch, creating new cache."); #endif - for(i = 0; i < currentCache->ringCount; i++) - { - freeTree(currentCache->ringIndices[i]); - } - lwfree(currentCache->ringIndices); - lwfree(currentCache->poly); - lwfree(currentCache); - return createNewCache(poly, serializedPoly); - } - for(i = 0; i < length; i++) - { - uchar a = serializedPoly[i]; - uchar b = currentCache->poly[i]; - if(a != b) - { + clearCache(currentCache); + return currentCache; + } + if( memcmp(serializedPoly, currentCache->poly, length) ) + { #ifdef PGIS_DEBUG - lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b); -#endif - for(i = 0; i < currentCache->ringCount; i++) - { - freeTree(currentCache->ringIndices[i]); - } - lwfree(currentCache->ringIndices); - lwfree(currentCache->poly); - lwfree(currentCache); - return createNewCache(poly, serializedPoly); - } - } + lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b); +#endif + clearCache(currentCache); + return currentCache; + } #ifdef PGIS_DEBUG lwnotice("Polygon match, retaining current cache, %p.", currentCache); diff --git a/lwgeom/lwgeom_rtree.h b/lwgeom/lwgeom_rtree.h index db17de7cc..cc694cfb2 100644 --- a/lwgeom/lwgeom_rtree.h +++ b/lwgeom/lwgeom_rtree.h @@ -47,6 +47,7 @@ typedef struct { RTREE_NODE **ringIndices; int ringCount; + int polyCount; uchar *poly; } RTREE_POLY_CACHE; @@ -54,6 +55,10 @@ typedef struct * Creates a new cachable index if needed, or returns the current cache if * it is applicable to the current polygon. */ -RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, RTREE_POLY_CACHE *currentCache); +RTREE_POLY_CACHE *retrieveCache(LWGEOM *lwgeom, uchar *serializedPoly, RTREE_POLY_CACHE *currentCache); +RTREE_POLY_CACHE *createCache(void); +/* Frees the cache. */ +void populateCache(RTREE_POLY_CACHE *cache, LWGEOM *lwgeom, uchar *serializedPoly); +void clearCache(RTREE_POLY_CACHE *cache); #endif /* !defined _LIBLWGEOM_H */ -- 2.40.0