From: Paul Ramsey Date: Fri, 26 Sep 2008 19:43:19 +0000 (+0000) Subject: P-I-P rennovation complete: memory leaks gone, multipolygon support added X-Git-Tag: 1.4.0b1~695 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=a4fa95b3017261c2940ea088c6b226ae75616773;p=postgis P-I-P rennovation complete: memory leaks gone, multipolygon support added git-svn-id: http://svn.osgeo.org/postgis/trunk@3017 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/lwgeom/lwgeom_functions_analytic.c b/lwgeom/lwgeom_functions_analytic.c index 48884719f..138054fed 100644 --- a/lwgeom/lwgeom_functions_analytic.c +++ b/lwgeom/lwgeom_functions_analytic.c @@ -1276,6 +1276,63 @@ int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) return 1; } +/* + * return -1 if point outside polygon + * 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. + */ +int point_in_multipolygon(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point) +{ + int i; + POINT2D pt; + int result = -1; + + LWDEBUGF(2, "point_in_multipolygon called for %p %d %d %p.", root, polyCount, ringCount, 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(root[i], &pt); + LWDEBUGF(4, "point_in_multipolygon: exterior ring (%d), point_in_ring returned %d", i, in_ring); + if( in_ring != -1 ) /* not outside this ring */ + { + LWDEBUG(3, "point_in_multipolygon: inside exterior ring."); + 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(root[i], &pt); + LWDEBUGF(4, "point_in_multipolygon: hole (%d), point_in_ring returned %d", i, in_ring); + if( in_ring == 1 ) /* completely inside hole */ + { + LWDEBUGF(3, "point_in_multipolygon: within hole %d.", i); + 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 outside polygon or on boundary * return 1 iff point inside polygon @@ -1349,6 +1406,7 @@ int point_outside_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point) return 0; } + /* * return 0 iff point inside polygon or on boundary * return 1 iff point outside polygon diff --git a/lwgeom/lwgeom_geos_c.c b/lwgeom/lwgeom_geos_c.c index 1a31f06b4..84d880ba5 100644 --- a/lwgeom/lwgeom_geos_c.c +++ b/lwgeom/lwgeom_geos_c.c @@ -1382,7 +1382,7 @@ 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); +int point_in_multipolygon(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point); PG_FUNCTION_INFO_V1(contains); @@ -1391,14 +1391,13 @@ Datum contains(PG_FUNCTION_ARGS) PG_LWGEOM *geom1; PG_LWGEOM *geom2; GEOSGeom g1,g2; - bool result; BOX2DFLOAT4 box1, box2; - int type1, type2; - LWPOLY *poly; - /* LWMPOLY *mpoly; */ - LWPOINT *point; - RTREE_POLY_CACHE *poly_cache; - MemoryContext old_context; + int type1, type2; + LWGEOM *lwgeom; + LWPOINT *point; + RTREE_POLY_CACHE *poly_cache; + MemoryContext old_context; + bool result; #ifdef PROFILE profstart(PROF_QRUN); @@ -1410,90 +1409,63 @@ Datum contains(PG_FUNCTION_ARGS) errorIfGeometryCollection(geom1,geom2); errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2)); - /* - * short-circuit 1: if geom2 bounding box is not completely inside - * 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) ) + /* + ** short-circuit 1: if geom2 bounding box is not completely inside + ** 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 ( 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 - * call the point-in-polygon function. - */ - type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); - type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); - if(type1 == POLYGONTYPE && type2 == POINTTYPE) - { - POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); - point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); + /* + ** short-circuit 2: if geom2 is a point and geom1 is a polygon + ** call the point-in-polygon function. + */ + type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]); + type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]); + if((type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE) && type2 == POINTTYPE) + { + POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); + point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); - POSTGIS_DEBUGF(3, "Precall point_in_polygon %p, %p", poly, point); + POSTGIS_DEBUGF(3, "Precall point_in_multipolygon %p, %p", lwgeom, point); - /* - * 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_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); - } - } - /* Not yet functional - else if(type1 == MULTIPOLYGONTYPE && type2 == POINTTYPE) - { - POSTGIS_DEBUG(3, "Point in MultiPolygon test requested...short-circuiting."); - - 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 - { - POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", type1, type2); - } + /* + * 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(geom1), fcinfo->flinfo->fn_extra); + fcinfo->flinfo->fn_extra = poly_cache; + MemoryContextSwitchTo(old_context); + + result = point_in_multipolygon(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + 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); + } + } + else + { + POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", type1, type2); + } initGEOS(lwnotice, lwnotice); @@ -1554,7 +1526,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; @@ -1590,14 +1562,14 @@ 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) { POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); - poly = lwpoly_deserialize(SERIALIZED_FORM(geom1)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); - POSTGIS_DEBUGF(3, "Precall point_in_polygon %p, %p", poly, point); + POSTGIS_DEBUGF(3, "Precall point_in_polygon %p, %p", lwgeom, point); /* * Switch the context to the function-scope context, @@ -1605,32 +1577,29 @@ 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); - } + result = point_in_multipolygon(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + 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 { POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", type1, type2); } - + initGEOS(lwnotice, lwnotice); #ifdef PROFILE @@ -1684,7 +1653,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; @@ -1719,12 +1688,12 @@ 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 || type1 == MULTIPOLYGONTYPE) && type1 == POINTTYPE) { POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); 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, @@ -1732,26 +1701,23 @@ 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); - } + result = point_in_multipolygon(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + 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); @@ -1812,7 +1778,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; @@ -1850,12 +1816,12 @@ 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) { POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); 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, @@ -1863,26 +1829,23 @@ 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); - } + result = point_in_multipolygon(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + 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); @@ -2011,18 +1974,18 @@ Datum crosses(PG_FUNCTION_ARGS) } - PG_FUNCTION_INFO_V1(intersects); Datum intersects(PG_FUNCTION_ARGS) { PG_LWGEOM *geom1; PG_LWGEOM *geom2; + uchar *serialized_poly; GEOSGeom g1,g2; bool result; BOX2DFLOAT4 box1, box2; int type1, type2; LWPOINT *point; - LWPOLY *poly; + LWGEOM *lwgeom; MemoryContext old_context; RTREE_POLY_CACHE *poly_cache; @@ -2056,77 +2019,44 @@ 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))) { POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); - 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, 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( type1 == POINTTYPE ) { + point = lwpoint_deserialize(SERIALIZED_FORM(geom1)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom2)); + serialized_poly = SERIALIZED_FORM(geom2); + } else { + point = lwpoint_deserialize(SERIALIZED_FORM(geom2)); + lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom1)); + serialized_poly = SERIALIZED_FORM(geom1); } - } - else if(type1 == POLYGONTYPE && type2 == POINTTYPE) - { - POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting."); - - 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) + /* + * 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_poly, fcinfo->flinfo->fn_extra); + fcinfo->flinfo->fn_extra = poly_cache; + MemoryContextSwitchTo(old_context); + + result = point_in_multipolygon(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point); + 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); } - } - initGEOS(lwnotice, lwnotice); #ifdef PROFILE @@ -2258,11 +2188,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); @@ -2288,80 +2213,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) - { - POSTGIS_DEBUG(3, "Point outside Polygon test requested...short-circuiting."); - - 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) - { - POSTGIS_DEBUG(3, "Point outside Polygon test requested...short-circuiting."); - - 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); @@ -3558,17 +3409,19 @@ Datum intersectsPrepared(PG_FUNCTION_ARGS); /* * GEOS prepared geometry is only available from GEOS 3.1 onwards */ -#if GEOS_VERNUM >= 31 +#if POSTGIS_GEOS_VERSION >= 31 #define PREPARED_GEOM 1 +#warning COMPILING PREPARED GEOMETRY #endif - #ifdef PREPARED_GEOM typedef struct { - int32 key; - GEOSPreparedGeometry * prepared_geom; + int32 key1; + int32 key2; + int32 argnum; + GEOSPreparedGeometry *prepared_geom; } PREPARED_GEOM_CACHE; @@ -3591,46 +3444,60 @@ typedef struct * key into cache * clear prepared cache */ -PREPARED_GEOM_CACHE * -get_prepared_geometry_cache( - PREPARED_GEOM_CACHE * cache, - uchar * serialized_geom, - int32 key ) +PREPARED_GEOM_CACHE +*get_prepared_geometry_cache( + PREPARED_GEOM_CACHE *cache, + PG_LWGEOM *serialized_geom1, PG_LWGEOM *serialized_geom2, + int32 key1, int32 key2 ) { - GEOSGeom g; if ( !cache ) { - LWDEBUGF(3, "get_prepared_geometry_cache: creating cache: %x", cache); - - cache = lwalloc( sizeof( PREPARED_GEOM_CACHE )); - + cache = lwalloc( sizeof(PREPARED_GEOM_CACHE) ); cache->prepared_geom = 0; - cache->key = key; + cache->argnum = 0; + cache->key1 = 0; + cache->key2 = 0; + LWDEBUGF(3, "get_prepared_geometry_cache: creating cache: %x", cache); } - else if ( cache->key == key ) + else if ( cache->key1 == key1 ) { if ( !cache->prepared_geom ) { - LWDEBUGF(3, "get_prepared_geometry_cache: preparing obj"); - - g = POSTGIS2GEOS( serialized_geom); - cache->prepared_geom = GEOSPrepare( g); + GEOSGeom g = POSTGIS2GEOS( serialized_geom1 ); + cache->prepared_geom = GEOSPrepare( g ); + cache->argnum = 1; + LWDEBUG(3, "get_prepared_geometry_cache: preparing obj in argument 1"); } else { - LWDEBUGF(3, "get_prepared_geometry_cache: prepared obj in cache"); + LWDEBUG(3, "get_prepared_geometry_cache: prepared obj 1 in cache"); } } - else + else if ( key2 && cache->key2 == key2 ) + { + if ( !cache->prepared_geom ) + { + GEOSGeom g = POSTGIS2GEOS( serialized_geom2 ); + cache->prepared_geom = GEOSPrepare( g ); + cache->argnum = 2; + LWDEBUG(3, "get_prepared_geometry_cache: preparing obj in argument 2"); + } + else + { + LWDEBUG(3, "get_prepared_geometry_cache: prepared obj 2 in cache"); + } + } + else if ( cache->prepared_geom ) { LWDEBUG(3, "get_prepared_geometry_cache: obj NOT in cache"); - - GEOSPreparedGeom_destroy( cache->prepared_geom); - + GEOSPreparedGeom_destroy( cache->prepared_geom ); cache->prepared_geom = 0; - cache->key = key; + cache->argnum = 0; } + + cache->key1 = key1; + cache->key2 = key2; return cache; } @@ -3644,22 +3511,19 @@ Datum containsPrepared(PG_FUNCTION_ARGS) #ifndef PREPARED_GEOM elog(ERROR,"Not implemented in this version!"); PG_RETURN_NULL(); /* never get here */ - #else PG_LWGEOM * geom1; PG_LWGEOM * geom2; - GEOSGeom g1, g2; GEOSPreparedGeometry * pg; bool result; BOX2DFLOAT4 box1, box2; PREPARED_GEOM_CACHE * prep_cache; MemoryContext old_context; - int32 surrogate_key; + int32 key1, key2; geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - - surrogate_key = PG_GETARG_INT32(2); + key1 = PG_GETARG_INT32(2); errorIfGeometryCollection(geom1,geom2); errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2)); @@ -3672,39 +3536,34 @@ Datum containsPrepared(PG_FUNCTION_ARGS) 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); } old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); prep_cache = fcinfo->flinfo->fn_extra; - - prep_cache = get_prepared_geometry_cache( prep_cache, geom1, surrogate_key); - + prep_cache = get_prepared_geometry_cache( prep_cache, geom1, 0, key1, 0 ); fcinfo->flinfo->fn_extra = prep_cache; MemoryContextSwitchTo(old_context); initGEOS(lwnotice, lwnotice); - g2 = POSTGIS2GEOS(geom2); - - if ( !prep_cache || !prep_cache->prepared_geom ) + if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 ) { - g1 = POSTGIS2GEOS(geom1); - - result = GEOSContains( g1, g2); - - GEOSGeom_destroy(g1); + GEOSGeom g = POSTGIS2GEOS(geom2); + pg = prep_cache->prepared_geom; + result = GEOSPreparedContains( pg, g); + GEOSGeom_destroy(g); } else { - pg = prep_cache->prepared_geom; - - result = GEOSPreparedContains( pg, g2); - } - GEOSGeom_destroy(g2); + GEOSGeom g1 = POSTGIS2GEOS(geom1); + GEOSGeom g2 = POSTGIS2GEOS(geom2); + result = GEOSContains( g1, g2); + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + } if (result == 2) { @@ -3725,22 +3584,20 @@ Datum containsProperlyPrepared(PG_FUNCTION_ARGS) #ifndef PREPARED_GEOM elog(ERROR,"Not implemented in this version!"); PG_RETURN_NULL(); /* never get here */ - #else PG_LWGEOM * geom1; PG_LWGEOM * geom2; - GEOSGeom g1, g2; GEOSPreparedGeometry * pg; bool result; BOX2DFLOAT4 box1, box2; PREPARED_GEOM_CACHE * prep_cache; MemoryContext old_context; - int32 surrogate_key; + int32 key1, key2; geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - - surrogate_key = PG_GETARG_INT32(2); + key1 = PG_GETARG_INT32(1); + geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2)); + key2 = PG_GETARG_INT32(3); errorIfGeometryCollection(geom1,geom2); errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2)); @@ -3753,39 +3610,34 @@ Datum containsProperlyPrepared(PG_FUNCTION_ARGS) 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); } old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); prep_cache = fcinfo->flinfo->fn_extra; - - prep_cache = get_prepared_geometry_cache( prep_cache, geom1, surrogate_key); - + prep_cache = get_prepared_geometry_cache( prep_cache, geom1, 0, key1, 0 ); fcinfo->flinfo->fn_extra = prep_cache; - MemoryContextSwitchTo( old_context); + MemoryContextSwitchTo(old_context); initGEOS(lwnotice, lwnotice); - - g2 = POSTGIS2GEOS(geom2); - if ( !prep_cache || !prep_cache->prepared_geom ) + if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 ) { - g1 = POSTGIS2GEOS(geom1); - - result = GEOSRelatePattern( g1, g2, "T**FF*FF*"); - - GEOSGeom_destroy(g1); + GEOSGeom g = POSTGIS2GEOS(geom2); + pg = prep_cache->prepared_geom; + result = GEOSPreparedContainsProperly( pg, g); + GEOSGeom_destroy(g); } else { - pg = prep_cache->prepared_geom; - - result = GEOSPreparedContainsProperly( pg, g2); - } - GEOSGeom_destroy(g2); + GEOSGeom g1 = POSTGIS2GEOS(geom1); + GEOSGeom g2 = POSTGIS2GEOS(geom2); + result = GEOSRelatePattern( g1, g2, "T**FF*FF*" ); + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + } if (result == 2) { @@ -3806,22 +3658,20 @@ Datum coversPrepared(PG_FUNCTION_ARGS) #ifndef PREPARED_GEOM elog(ERROR,"Not implemented in this version!"); PG_RETURN_NULL(); /* never get here */ - #else PG_LWGEOM * geom1; PG_LWGEOM * geom2; - GEOSGeom g1, g2; GEOSPreparedGeometry * pg; bool result; BOX2DFLOAT4 box1, box2; PREPARED_GEOM_CACHE * prep_cache; MemoryContext old_context; - int32 surrogate_key; + int32 key1, key2; geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - - surrogate_key = PG_GETARG_INT32(2); + key1 = PG_GETARG_INT32(1); + geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2)); + key2 = PG_GETARG_INT32(3); errorIfGeometryCollection(geom1,geom2); errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2)); @@ -3834,39 +3684,34 @@ Datum coversPrepared(PG_FUNCTION_ARGS) 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); } old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); prep_cache = fcinfo->flinfo->fn_extra; - - prep_cache = get_prepared_geometry_cache( prep_cache, geom1, surrogate_key); - + prep_cache = get_prepared_geometry_cache( prep_cache, geom1, 0, key1, 0 ); fcinfo->flinfo->fn_extra = prep_cache; MemoryContextSwitchTo(old_context); initGEOS(lwnotice, lwnotice); - g2 = POSTGIS2GEOS(geom2); - - if ( !prep_cache || !prep_cache->prepared_geom ) + if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 ) { - g1 = POSTGIS2GEOS(geom1); - - result = GEOSRelatePattern( g1, g2, "******FF*"); - - GEOSGeom_destroy(g1); + GEOSGeom g = POSTGIS2GEOS(geom2); + pg = prep_cache->prepared_geom; + result = GEOSPreparedCovers( pg, g); + GEOSGeom_destroy(g); } else { - pg = prep_cache->prepared_geom; - - result = GEOSPreparedCovers( pg, g2); - } - GEOSGeom_destroy(g2); + GEOSGeom g1 = POSTGIS2GEOS(geom1); + GEOSGeom g2 = POSTGIS2GEOS(geom2); + result = GEOSRelatePattern( g1, g2, "******FF*" ); + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + } if (result == 2) { @@ -3881,28 +3726,27 @@ Datum coversPrepared(PG_FUNCTION_ARGS) #endif /* PREPARED_GEOM */ } + PG_FUNCTION_INFO_V1(intersectsPrepared); Datum intersectsPrepared(PG_FUNCTION_ARGS) { #ifndef PREPARED_GEOM elog(ERROR,"Not implemented in this version!"); PG_RETURN_NULL(); /* never get here */ - #else PG_LWGEOM * geom1; PG_LWGEOM * geom2; - GEOSGeom g1, g2; GEOSPreparedGeometry * pg; bool result; BOX2DFLOAT4 box1, box2; PREPARED_GEOM_CACHE * prep_cache; MemoryContext old_context; - int32 surrogate_key; + int32 key1, key2; geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - - surrogate_key = PG_GETARG_INT32(2); + key1 = PG_GETARG_INT32(1); + geom2 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2)); + key2 = PG_GETARG_INT32(3); errorIfGeometryCollection(geom1,geom2); errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2)); @@ -3915,39 +3759,42 @@ Datum intersectsPrepared(PG_FUNCTION_ARGS) if ( getbox2d_p(SERIALIZED_FORM(geom1), &box1) && getbox2d_p(SERIALIZED_FORM(geom2), &box2) ) { - if ( box2.xmax < box1.xmin ) PG_RETURN_BOOL(FALSE); - if ( box2.xmin > box1.xmax ) PG_RETURN_BOOL(FALSE); - if ( box2.ymax < box1.ymin ) PG_RETURN_BOOL(FALSE); - if ( box2.ymin > box2.ymax ) PG_RETURN_BOOL(FALSE); + if (( box2.xmax < box1.xmin ) || ( box2.xmin > box1.xmax ) || + ( box2.ymax < box1.ymin ) || ( box2.ymin > box2.ymax )) + PG_RETURN_BOOL(FALSE); } old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt); prep_cache = fcinfo->flinfo->fn_extra; - - prep_cache = get_prepared_geometry_cache( prep_cache, geom1, surrogate_key); - + prep_cache = get_prepared_geometry_cache( prep_cache, geom1, geom2, key1, key2 ); fcinfo->flinfo->fn_extra = prep_cache; MemoryContextSwitchTo(old_context); initGEOS(lwnotice, lwnotice); - g2 = POSTGIS2GEOS(geom2); - - if ( !prep_cache || !prep_cache->prepared_geom ) + if ( prep_cache && prep_cache->prepared_geom ) { - g1 = POSTGIS2GEOS(geom1); - - result = GEOSIntersects( g1, g2); - - GEOSGeom_destroy(g1); + if( prep_cache->argnum == 1 ) { + GEOSGeom g = POSTGIS2GEOS(geom2); + pg = prep_cache->prepared_geom; + result = GEOSPreparedCovers( pg, g); + GEOSGeom_destroy(g); + } + else { + GEOSGeom g = POSTGIS2GEOS(geom1); + pg = prep_cache->prepared_geom; + result = GEOSPreparedCovers( pg, g); + GEOSGeom_destroy(g); + } } else { - pg = prep_cache->prepared_geom; - - result = GEOSPreparedIntersects( pg, g2); + GEOSGeom g1 = POSTGIS2GEOS(geom1); + GEOSGeom g2 = POSTGIS2GEOS(geom2); + result = GEOSIntersects( g1, g2); + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); } - GEOSGeom_destroy(g2); if (result == 2) { diff --git a/lwgeom/lwgeom_rtree.c b/lwgeom/lwgeom_rtree.c index dab5306be..b6aaa14da 100644 --- a/lwgeom/lwgeom_rtree.c +++ b/lwgeom/lwgeom_rtree.c @@ -24,80 +24,80 @@ Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS); */ RTREE_NODE *createTree(POINTARRAY *pointArray) { - RTREE_NODE *root; + RTREE_NODE *root; RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints); int i, nodeCount; - int childNodes, parentNodes; - - LWDEBUGF(2, "createTree called with pointarray %p", pointArray); - - nodeCount = pointArray->npoints - 1; - - LWDEBUGF(3, "Total leaf nodes: %d", nodeCount); - - /* - * Create a leaf node for every line segment. - */ - for(i = 0; i < nodeCount; i++) - { - nodes[i] = createLeafNode(pointArray, i); - } - - /* - * Next we group nodes by pairs. If there's an odd number of nodes, - * we bring the last node up a level as is. Continue until we have - * a single top node. - */ - childNodes = nodeCount; - parentNodes = nodeCount / 2; - while(parentNodes > 0) - { - LWDEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes); - - i = 0; - while(i < parentNodes) - { - nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]); - i++; - } - /* - * Check for an odd numbered final node. - */ - if(parentNodes * 2 < childNodes) - { - LWDEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i); - - nodes[i] = nodes[childNodes - 1]; - parentNodes++; - } - childNodes = parentNodes; - parentNodes = parentNodes / 2; - } - - root = nodes[0]; - lwfree(nodes); - LWDEBUGF(3, "createTree returning %p", root); - - return root; + int childNodes, parentNodes; + + LWDEBUGF(2, "createTree called with pointarray %p", pointArray); + + nodeCount = pointArray->npoints - 1; + + LWDEBUGF(3, "Total leaf nodes: %d", nodeCount); + + /* + * Create a leaf node for every line segment. + */ + for(i = 0; i < nodeCount; i++) + { + nodes[i] = createLeafNode(pointArray, i); + } + + /* + * Next we group nodes by pairs. If there's an odd number of nodes, + * we bring the last node up a level as is. Continue until we have + * a single top node. + */ + childNodes = nodeCount; + parentNodes = nodeCount / 2; + while(parentNodes > 0) + { + LWDEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes); + + i = 0; + while(i < parentNodes) + { + nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]); + i++; + } + /* + * Check for an odd numbered final node. + */ + if(parentNodes * 2 < childNodes) + { + LWDEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i); + + nodes[i] = nodes[childNodes - 1]; + parentNodes++; + } + childNodes = parentNodes; + parentNodes = parentNodes / 2; + } + + root = nodes[0]; + lwfree(nodes); + LWDEBUGF(3, "createTree returning %p", root); + + return root; } /* * Creates an interior node given the children. */ RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){ - RTREE_NODE *parent; + RTREE_NODE *parent; - LWDEBUGF(2, "createInteriorNode called for children %p, %p", left, right); + LWDEBUGF(2, "createInteriorNode called for children %p, %p", left, right); - parent = lwalloc(sizeof(RTREE_NODE)); - parent->leftNode = left; - parent->rightNode = right; - parent->interval = mergeIntervals(left->interval, right->interval); - parent->segment = NULL; + parent = lwalloc(sizeof(RTREE_NODE)); + parent->leftNode = left; + parent->rightNode = right; + parent->interval = mergeIntervals(left->interval, right->interval); + parent->segment = NULL; - LWDEBUGF(3, "createInteriorNode returning %p", parent); + LWDEBUGF(3, "createInteriorNode returning %p", parent); - return parent; + return parent; } /* @@ -105,52 +105,52 @@ RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){ */ RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint) { - RTREE_NODE *parent; - LWLINE *line; - double value1; - double value2; - POINT4D tmp; - POINTARRAY *npa; - - LWDEBUGF(2, "createLeafNode called for point %d of %p", startPoint, pa); - - if(pa->npoints < startPoint + 2) - { - lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint); - } - - /* - * The given point array will be part of a geometry that will be freed - * independently of the index. Since we may want to cache the index, - * we must create independent arrays. - */ - npa = lwalloc(sizeof(POINTARRAY)); - npa->dims = 0; - npa->npoints = 2; - TYPE_SETZM(npa->dims, 0, 0); - npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2); - - getPoint4d_p(pa, startPoint, &tmp); - - setPoint4d(npa, 0, &tmp); - value1 = tmp.y; - - getPoint4d_p(pa, startPoint + 1, &tmp); - - setPoint4d(npa, 1, &tmp); - value2 = tmp.y; - - line = lwline_construct(-1, NULL, npa); - - parent = lwalloc(sizeof(RTREE_NODE)); - parent->interval = createInterval(value1, value2); - parent->segment = line; - parent->leftNode = NULL; - parent->rightNode = NULL; - - LWDEBUGF(3, "createLeafNode returning %p", parent); - - return parent; + RTREE_NODE *parent; + LWLINE *line; + double value1; + double value2; + POINT4D tmp; + POINTARRAY *npa; + + LWDEBUGF(2, "createLeafNode called for point %d of %p", startPoint, pa); + + if(pa->npoints < startPoint + 2) + { + lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint); + } + + /* + * The given point array will be part of a geometry that will be freed + * independently of the index. Since we may want to cache the index, + * we must create independent arrays. + */ + npa = lwalloc(sizeof(POINTARRAY)); + npa->dims = 0; + npa->npoints = 2; + TYPE_SETZM(npa->dims, 0, 0); + npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2); + + getPoint4d_p(pa, startPoint, &tmp); + + setPoint4d(npa, 0, &tmp); + value1 = tmp.y; + + getPoint4d_p(pa, startPoint + 1, &tmp); + + setPoint4d(npa, 1, &tmp); + value2 = tmp.y; + + line = lwline_construct(-1, NULL, npa); + + parent = lwalloc(sizeof(RTREE_NODE)); + parent->interval = createInterval(value1, value2); + parent->segment = line; + parent->leftNode = NULL; + parent->rightNode = NULL; + + LWDEBUGF(3, "createLeafNode returning %p", parent); + + return parent; } /* @@ -158,17 +158,17 @@ RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint) */ INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2) { - INTERVAL *interval; + INTERVAL *interval; - LWDEBUGF(2, "mergeIntervals called with %p, %p", inter1, inter2); + LWDEBUGF(2, "mergeIntervals called with %p, %p", inter1, inter2); - interval = lwalloc(sizeof(INTERVAL)); - interval->max = FP_MAX(inter1->max, inter2->max); - interval->min = FP_MIN(inter1->min, inter2->min); + interval = lwalloc(sizeof(INTERVAL)); + interval->max = FP_MAX(inter1->max, inter2->max); + interval->min = FP_MIN(inter1->min, inter2->min); - LWDEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max); + LWDEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max); - return interval; + return interval; } /* @@ -176,17 +176,17 @@ INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2) */ INTERVAL *createInterval(double value1, double value2) { - INTERVAL *interval; + INTERVAL *interval; - LWDEBUGF(2, "createInterval called with %8.3f, %8.3f", value1, value2); + LWDEBUGF(2, "createInterval called with %8.3f, %8.3f", value1, value2); - interval = lwalloc(sizeof(INTERVAL)); - interval->max = FP_MAX(value1, value2); - interval->min = FP_MIN(value1, value2); + interval = lwalloc(sizeof(INTERVAL)); + interval->max = FP_MAX(value1, value2); + interval->min = FP_MIN(value1, value2); - LWDEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max); + LWDEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max); - return interval; + return interval; } /* @@ -195,19 +195,19 @@ INTERVAL *createInterval(double value1, double value2) */ void freeTree(RTREE_NODE *root) { - LWDEBUGF(2, "freeTree called for %p", 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); + LWDEBUGF(2, "freeTree called for %p", 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); } @@ -216,15 +216,15 @@ void freeTree(RTREE_NODE *root) */ void freeCache(RTREE_POLY_CACHE *cache) { - int i; - LWDEBUGF(2, "freeCache called for %p", cache); - for(i = 0; i < cache->ringCount; i++) - { - freeTree(cache->ringIndices[i]); - } - lwfree(cache->ringIndices); - lwfree(cache->poly); - lwfree(cache); + int i; + LWDEBUGF(2, "freeCache called for %p", cache); + for(i = 0; i < cache->ringCount; i++) + { + freeTree(cache->ringIndices[i]); + } + lwfree(cache->ringIndices); + lwfree(cache->poly); + lwfree(cache); } @@ -236,222 +236,275 @@ void freeCache(RTREE_POLY_CACHE *cache) */ LWMLINE *findLineSegments(RTREE_NODE *root, double value) { - LWMLINE *tmp, *result; - LWGEOM **lwgeoms; - - LWDEBUGF(2, "findLineSegments called for tree %p and value %8.3f", root, value); - - result = NULL; - - if(!isContained(root->interval, value)) - { - LWDEBUGF(3, "findLineSegments %p: not contained.", root); - - return NULL; - } - - /* If there is a segment defined for this node, include it. */ - if(root->segment) - { - LWDEBUGF(3, "findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type)); - - lwgeoms = lwalloc(sizeof(LWGEOM *)); - lwgeoms[0] = (LWGEOM *)root->segment; - - LWDEBUGF(3, "Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type)); - - result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms); - } - - /* If there is a left child node, recursively include its results. */ - if(root->leftNode) - { - LWDEBUGF(3, "findLineSegments %p: recursing left.", root); - - tmp = findLineSegments(root->leftNode, value); - if(tmp) - { - LWDEBUGF(3, "Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); - - if(result) - result = mergeMultiLines(result, tmp); - else - result = tmp; - } - } - - /* Same for any right child. */ - if(root->rightNode) - { - LWDEBUGF(3, "findLineSegments %p: recursing right.", root); - - tmp = findLineSegments(root->rightNode, value); - if(tmp) - { - LWDEBUGF(3, "Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); - - if(result) - result = mergeMultiLines(result, tmp); - else - result = tmp; - } - } - - return result; + LWMLINE *tmp, *result; + LWGEOM **lwgeoms; + + LWDEBUGF(2, "findLineSegments called for tree %p and value %8.3f", root, value); + + result = NULL; + + if(!isContained(root->interval, value)) + { + LWDEBUGF(3, "findLineSegments %p: not contained.", root); + + return NULL; + } + + /* If there is a segment defined for this node, include it. */ + if(root->segment) + { + LWDEBUGF(3, "findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type)); + + lwgeoms = lwalloc(sizeof(LWGEOM *)); + lwgeoms[0] = (LWGEOM *)root->segment; + + LWDEBUGF(3, "Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type)); + + result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms); + } + + /* If there is a left child node, recursively include its results. */ + if(root->leftNode) + { + LWDEBUGF(3, "findLineSegments %p: recursing left.", root); + + tmp = findLineSegments(root->leftNode, value); + if(tmp) + { + LWDEBUGF(3, "Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); + + if(result) + result = mergeMultiLines(result, tmp); + else + result = tmp; + } + } + + /* Same for any right child. */ + if(root->rightNode) + { + LWDEBUGF(3, "findLineSegments %p: recursing right.", root); + + tmp = findLineSegments(root->rightNode, value); + if(tmp) + { + LWDEBUGF(3, "Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); + + if(result) + result = mergeMultiLines(result, tmp); + else + result = tmp; + } + } + + return result; } /* Merges two multilinestrings into a single multilinestring. */ LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2) { - LWGEOM **geoms; - LWCOLLECTION *col; - int i, j, ngeoms; - - LWDEBUGF(2, "mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type)); - - ngeoms = line1->ngeoms + line2->ngeoms; - geoms = lwalloc(sizeof(LWGEOM *) * ngeoms); - - j = 0; - for(i = 0; i < line1->ngeoms; i++, j++) - { - geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]); - } - for(i = 0; i < line2->ngeoms; i++, j++) - { - geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]); - } - col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms); - - LWDEBUGF(3, "mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type)); - - return (LWMLINE *)col; + LWGEOM **geoms; + LWCOLLECTION *col; + int i, j, ngeoms; + + LWDEBUGF(2, "mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type)); + + ngeoms = line1->ngeoms + line2->ngeoms; + geoms = lwalloc(sizeof(LWGEOM *) * ngeoms); + + j = 0; + for(i = 0; i < line1->ngeoms; i++, j++) + { + geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]); + } + for(i = 0; i < line2->ngeoms; i++, j++) + { + geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]); + } + col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms); + + LWDEBUGF(3, "mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type)); + + return (LWMLINE *)col; } /* * Returns 1 if min < value <= max, 0 otherwise. */ uint32 isContained(INTERVAL *interval, double value) { - return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0; + return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0; } PG_FUNCTION_INFO_V1(LWGEOM_polygon_index); Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS) { - PG_LWGEOM *igeom, *result; - LWGEOM *geom; - LWPOLY *poly; - LWMLINE *mline; - RTREE_NODE *root; - double yval; + PG_LWGEOM *igeom, *result; + LWGEOM *geom; + LWPOLY *poly; + LWMLINE *mline; + RTREE_NODE *root; + double yval; #if POSTGIS_DEBUG_LEVEL > 0 int i = 0; -#endif +#endif - POSTGIS_DEBUG(2, "polygon_index called."); - - result = NULL; - igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - yval = PG_GETARG_FLOAT8(1); - geom = lwgeom_deserialize(SERIALIZED_FORM(igeom)); - if(TYPE_GETTYPE(geom->type) != POLYGONTYPE) - { - lwgeom_release(geom); - PG_FREE_IF_COPY(igeom, 0); - PG_RETURN_NULL(); - } - poly = (LWPOLY *)geom; - root = createTree(poly->rings[0]); - - mline = findLineSegments(root, yval); + POSTGIS_DEBUG(2, "polygon_index called."); + + result = NULL; + igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + yval = PG_GETARG_FLOAT8(1); + geom = lwgeom_deserialize(SERIALIZED_FORM(igeom)); + if(TYPE_GETTYPE(geom->type) != POLYGONTYPE) + { + lwgeom_release(geom); + PG_FREE_IF_COPY(igeom, 0); + PG_RETURN_NULL(); + } + poly = (LWPOLY *)geom; + root = createTree(poly->rings[0]); + + mline = findLineSegments(root, yval); #if POSTGIS_DEBUG_LEVEL >= 3 - POSTGIS_DEBUGF(3, "mline returned %p %d", mline, TYPE_GETTYPE(mline->type)); - for(i = 0; i < mline->ngeoms; i++) - { - POSTGIS_DEBUGF(3, "geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type)); - } + POSTGIS_DEBUGF(3, "mline returned %p %d", mline, TYPE_GETTYPE(mline->type)); + for(i = 0; i < mline->ngeoms; i++) + { + POSTGIS_DEBUGF(3, "geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type)); + } #endif - if(mline) - result = pglwgeom_serialize((LWGEOM *)mline); + if(mline) + result = pglwgeom_serialize((LWGEOM *)mline); - POSTGIS_DEBUGF(3, "returning result %p", result); + POSTGIS_DEBUGF(3, "returning result %p", result); lwfree(root); PG_FREE_IF_COPY(igeom, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)mline); - PG_RETURN_POINTER(result); - + lwgeom_release((LWGEOM *)poly); + lwgeom_release((LWGEOM *)mline); + PG_RETURN_POINTER(result); + } -RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly) -{ - RTREE_POLY_CACHE *result; - int i, length; - - LWDEBUGF(2, "createNewCache called with %p", poly); - - 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]); - } - LWDEBUGF(3, "createNewCache returning %p", result); - - return result; +RTREE_POLY_CACHE *createNewCache(LWGEOM *lwgeom, uchar *serializedPoly) +{ + RTREE_POLY_CACHE *result; + int i, j, k, length; + + LWDEBUGF(2, "createNewCache called with %p", lwgeom); + + result = lwalloc(sizeof(RTREE_POLY_CACHE)); + + if(TYPE_GETTYPE(lwgeom->type) == MULTIPOLYGONTYPE) + { + LWMPOLY *mpoly = (LWMPOLY *)lwgeom; + int nrings = 0; + /* + ** Count the total number of rings. + */ + for( i = 0; i < mpoly->ngeoms; i++ ) + { + nrings += mpoly->geoms[i]->nrings; + } + result->polyCount = mpoly->ngeoms; + result->ringCount = nrings; + result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings); + /* + ** Load the exterior rings onto the ringIndices array first + */ + for( i = 0; i < mpoly->ngeoms; i++ ) + { + result->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++ ) + { + result->ringIndices[i] = createTree(mpoly->geoms[j]->rings[k]); + i++; + } + } + } + else if ( TYPE_GETTYPE(lwgeom->type) == POLYGONTYPE ) + { + LWPOLY *poly = (LWPOLY *)lwgeom; + result->polyCount = 1; + result->ringCount = poly->nrings; + /* + ** Just load the rings on in order + */ + result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings); + for( i = 0; i < poly->nrings; i++ ) + { + result->ringIndices[i] = createTree(poly->rings[i]); + } + } + else + { + /* Uh oh, shouldn't be here. */ + return NULL; + } + + /* + ** Copy the serialized form of the polygon into the cache so + ** we can test for equality against subsequent polygons. + */ + length = lwgeom_size(serializedPoly); + result->poly = lwalloc(length); + memcpy(result->poly, serializedPoly, length); + + LWDEBUGF(3, "createNewCache returning %p", result); + + return result; } /* * Creates a new cachable index if needed, or returns the current cache if * it is applicable to the current polygon. * The memory context must be changed to function scope before calling this - * method. The method will allocate memory for the cache it creates, + * 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 *currentCache) +RTREE_POLY_CACHE *retrieveCache(LWGEOM *lwgeom, uchar *serializedPoly, + RTREE_POLY_CACHE *currentCache) { - int length; - - LWDEBUGF(2, "retrieveCache called with %p %p %p", poly, serializedPoly, currentCache); - - if(!currentCache) - { - LWDEBUG(3, "No existing cache, create one."); - return createNewCache(poly, serializedPoly); - } - if(!(currentCache->poly)) - { - LWDEBUG(3, "Cache contains no polygon, creating new cache."); - return createNewCache(poly, serializedPoly); - } - - length = lwgeom_size_poly(serializedPoly); - - if(lwgeom_size_poly(currentCache->poly) != length) - { - LWDEBUG(3, "Polygon size mismatch, creating new cache."); - freeCache(currentCache); - return createNewCache(poly, serializedPoly); - } - if( memcmp(serializedPoly, currentCache->poly, length) ) - { - LWDEBUG(3, "Polygon mismatch, creating new cache."); - freeCache(currentCache); - return createNewCache(poly, serializedPoly); - } - - LWDEBUGF(3, "Polygon match, retaining current cache, %p.", currentCache); - - return currentCache; + int length; + + LWDEBUGF(2, "retrieveCache called with %p %p %p", lwgeom, serializedPoly, currentCache); + + if(!currentCache) + { + LWDEBUG(3, "No existing cache, create one."); + return createNewCache(lwgeom, serializedPoly); + } + if(!(currentCache->poly)) + { + LWDEBUG(3, "Cache contains no polygon, creating new cache."); + return createNewCache(lwgeom, serializedPoly); + } + + length = lwgeom_size(serializedPoly); + + if(lwgeom_size(currentCache->poly) != length) + { + LWDEBUG(3, "Polygon size mismatch, creating new cache."); + freeCache(currentCache); + return createNewCache(lwgeom, serializedPoly); + } + if( memcmp(serializedPoly, currentCache->poly, length) ) + { + LWDEBUG(3, "Polygon mismatch, creating new cache."); + freeCache(currentCache); + return createNewCache(lwgeom, serializedPoly); + } + + LWDEBUGF(3, "Polygon match, retaining current cache, %p.", currentCache); + + return currentCache; } diff --git a/lwgeom/lwgeom_rtree.h b/lwgeom/lwgeom_rtree.h index aa82c3ee6..9413a9d08 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,8 +55,8 @@ 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 *createNewCache(LWPOLY *poly, uchar *serializedPoly); +RTREE_POLY_CACHE *retrieveCache(LWGEOM *lwgeom, uchar *serializedPoly, RTREE_POLY_CACHE *currentCache); +RTREE_POLY_CACHE *createNewCache(LWGEOM *lwgeom, uchar *serializedPoly); /* Frees the cache. */ void freeCache(RTREE_POLY_CACHE *cache);