#include "math.h"
#include "lwgeom_rtree.h"
#include "lwalgorithm.h"
+#include "lwgeom_functions_analytic.h"
/***********************************************************************
double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
int point_in_ring(POINTARRAY *pts, POINT2D *point);
-int point_in_polygon(LWPOLY *polygon, LWPOINT *point);
-int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point);
int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point);
-int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point);
-int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point);
PG_FUNCTION_INFO_V1(LWGEOM_simplify2d);
* return 0 if point on boundary
* return 1 if point inside polygon
*
- * Expected **root order is all the exterior rings first, then all the holes
- *
- * TODO: this could be made slightly more efficient by ordering the rings in
- * EIIIEIIIEIEI order (exterior/interior) and including list of exterior ring
- * positions on the cache object.
+ * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI
*/
-int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point)
+int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
{
- int i;
+ int i, p, r, in_ring;
POINT2D pt;
int result = -1;
- LWDEBUGF(2, "point_in_multipolygon_rtree called for %p %d %d %p.", root, polyCount, ringCount, point);
+ LWDEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point);
getPoint2d_p(point->point, 0, &pt);
/* assume bbox short-circuit has already been attempted */
- /* is the point inside (not outside) any of the exterior rings? */
- for ( i = 0; i < polyCount; i++ )
- {
- int in_ring = point_in_ring_rtree(root[i], &pt);
- LWDEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", i, in_ring);
- if ( in_ring != -1 ) /* not outside this ring */
- {
- LWDEBUG(3, "point_in_multipolygon_rtree: inside exterior ring.");
- result = in_ring;
- break;
- }
- }
-
- if ( result == -1 ) /* strictly outside all rings */
- return result;
+ i = 0; /* the current index into the root array */
- /* ok, it's in a ring, but if it's in a hole it's still outside */
- for ( i = polyCount; i < ringCount; i++ )
+ /* is the point inside any of the sub-polygons? */
+ for ( p = 0; p < polyCount; p++ )
{
- int in_ring = point_in_ring_rtree(root[i], &pt);
- LWDEBUGF(4, "point_in_multipolygon_rtree: hole (%d), point_in_ring returned %d", i, in_ring);
- if ( in_ring == 1 ) /* completely inside hole */
- {
- LWDEBUGF(3, "point_in_multipolygon_rtree: within hole %d.", i);
- return -1;
- }
- if ( in_ring == 0 ) /* on the boundary of a hole */
+ in_ring = point_in_ring_rtree(root[i], &pt);
+ LWDEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring);
+ if ( in_ring == -1 ) /* outside the exterior ring */
{
- result = 0;
+ LWDEBUG(3, "point_in_multipolygon_rtree: outside exterior ring.");
+ continue;
}
+ if ( in_ring == 0 ) /* on the boundary */
+ {
+ LWDEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p);
+ return 0;
+ }
+
+ result = in_ring;
+
+ for(r=1; r<ringCounts[p]; r++)
+ {
+ in_ring = point_in_ring_rtree(root[i+r], &pt);
+ LWDEBUGF(4, "point_in_multipolygon_rtree: interior ring (%d), point_in_ring returned %d", r, in_ring);
+ if (in_ring == 1) /* inside a hole => outside the polygon */
+ {
+ LWDEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p);
+ result = -1;
+ break;
+ }
+ if (in_ring == 0) /* on the edge of a hole */
+ {
+ LWDEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p);
+ return 0;
+ }
+ }
+ /* if we have a positive result, we can short-circuit and return it */
+ if ( result != -1)
+ {
+ return result;
+ }
+ /* increment the index by the total number of rings in the sub-poly */
+ /* we do this here in case we short-cutted out of the poly before looking at all the rings */
+ i += ringCounts[p];
}
+
return result; /* -1 = outside, 0 = boundary, 1 = inside */
}
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2001-2011 Refractions Research Inc.
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "lwgeom_rtree.h"
+
+/*
+** Public prototypes for analytic functions.
+*/
+
+int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point);
+int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point);
+int point_in_polygon(LWPOLY *polygon, LWPOINT *point);
+int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *pont);
+
#include "lwgeom_rtree.h"
#include "lwgeom_geos_prepared.h"
#include "funcapi.h"
+#include "lwgeom_functions_analytic.h"
#include <string.h>
#include <assert.h>
Datum pgis_union_geometry_array(PG_FUNCTION_ARGS);
-/** @todo TODO: move these to a lwgeom_functions_analytic.h
- */
-int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point);
-int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point);
-int point_in_polygon(LWPOLY *polygon, LWPOINT *point);
-int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *pont);
-
/*
** Prototypes end
*/
if ( poly_cache->ringIndices )
{
- result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point);
+ result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
}
else if ( type1 == POLYGONTYPE )
{
if ( poly_cache->ringIndices )
{
- result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point);
+ result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
}
else if ( type1 == POLYGONTYPE )
{
if ( poly_cache->ringIndices )
{
- result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point);
+ result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
}
else if ( type2 == POLYGONTYPE )
{
if ( poly_cache->ringIndices )
{
- result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point);
+ result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
}
else if ( type2 == POLYGONTYPE )
{
if ( poly_cache->ringIndices )
{
- result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCount, point);
+ result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
}
else if ( polytype == POLYGONTYPE )
{
*/
void clearCache(RTREE_POLY_CACHE *cache)
{
- int i;
+ int g, r, i;
LWDEBUGF(2, "clearCache called for %p", cache);
- for (i = 0; i < cache->ringCount; i++)
- {
- freeTree(cache->ringIndices[i]);
+ i = 0;
+ for (g = 0; g < cache->polyCount; g++)
+ {
+ for (r = 0; r < cache->ringCounts; r++)
+ {
+ freeTree(cache->ringIndices[i]);
+ i++;
+ }
}
lwfree(cache->ringIndices);
+ lwfree(cache->ringCounts);
lwfree(cache->poly);
cache->poly = 0;
cache->ringIndices = 0;
- cache->ringCount = 0;
+ cache->ringCounts = 0;
cache->polyCount = 0;
}
RTREE_POLY_CACHE *result;
result = lwalloc(sizeof(RTREE_POLY_CACHE));
result->polyCount = 0;
- result->ringCount = 0;
+ result->ringCounts = 0;
result->ringIndices = 0;
result->poly = 0;
result->type = 1;
void populateCache(RTREE_POLY_CACHE *currentCache, LWGEOM *lwgeom, PG_LWGEOM *serializedPoly)
{
- int i, j, k, length;
+ int i, p, r, length;
LWMPOLY *mpoly;
LWPOLY *poly;
int nrings;
/*
** Count the total number of rings.
*/
+ currentCache->polyCount = mpoly->ngeoms;
+ currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms);
for ( i = 0; i < mpoly->ngeoms; i++ )
{
+ currentCache->ringCounts[i] = mpoly->geoms[i]->nrings;
nrings += mpoly->geoms[i]->nrings;
}
- currentCache->polyCount = mpoly->ngeoms;
- currentCache->ringCount = nrings;
currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings);
/*
- ** Load the exterior rings onto the ringIndices array first
- */
- for ( i = 0; i < mpoly->ngeoms; i++ )
- {
- currentCache->ringIndices[i] = createTree(mpoly->geoms[i]->rings[0]);
- }
- /*
- ** Load the interior rings (holes) onto ringIndices next
+ ** Load the array in geometry order, each outer ring followed by the inner rings
+ ** associated with that outer ring
*/
- for ( j = 0; j < mpoly->ngeoms; j++ )
+ i = 0;
+ for ( p = 0; p < mpoly->ngeoms; p++ )
{
- for ( k = 1; k < mpoly->geoms[j]->nrings; k++ )
+ for ( r = 0; r < mpoly->geoms[p]->nrings; r++ )
{
- currentCache->ringIndices[i] = createTree(mpoly->geoms[j]->rings[k]);
+ currentCache->ringIndices[i] = createTree(mpoly->geoms[p]->rings[r]);
i++;
}
}
LWDEBUG(2, "populateCache POLYGON");
poly = (LWPOLY *)lwgeom;
currentCache->polyCount = 1;
- currentCache->ringCount = poly->nrings;
+ currentCache->ringCounts = lwalloc(sizeof(int));
+ currentCache->ringCounts[0] = poly->nrings;
/*
** Just load the rings on in order
*/
{\r
char type;\r
RTREE_NODE **ringIndices;\r
- int ringCount;\r
+ int* ringCounts;\r
int polyCount;\r
PG_LWGEOM *poly;\r
}\r
projects / postgis / commitdiff