lwgeom_rtree.c 12.6 KB
Newer Older
1 2 3
/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
4
 * http://postgis.net
5
 *
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright (C) 2001-2005 Refractions Research Inc.
22 23 24
 *
 **********************************************************************/

25

26 27 28
#include <assert.h>

#include "../postgis_config.h"
29 30
#include "lwgeom_pg.h"
#include "liblwgeom.h"
31
#include "liblwgeom_internal.h"         /* For FP comparators. */
32
#include "lwgeom_cache.h"
33 34 35
#include "lwgeom_rtree.h"


36 37
/* Prototypes */
static void RTreeFree(RTREE_NODE* root);
38 39

/**
40 41
* Allocate a fresh clean RTREE_POLY_CACHE
*/
42
static RTREE_POLY_CACHE*
43
RTreeCacheCreate()
44
{
45 46 47 48 49
	RTREE_POLY_CACHE* result;
	result = lwalloc(sizeof(RTREE_POLY_CACHE));
	memset(result, 0, sizeof(RTREE_POLY_CACHE));
	return result;
}
50

51 52 53 54
/**
* Recursively frees the child nodes, the interval and the line before
* freeing the root node.
*/
55
static void
56 57 58
RTreeFree(RTREE_NODE* root)
{
	POSTGIS_DEBUGF(2, "RTreeFree called for %p", root);
59

60 61 62 63 64 65
	if (root->leftNode)
		RTreeFree(root->leftNode);
	if (root->rightNode)
		RTreeFree(root->rightNode);
	lwfree(root->interval);
	if (root->segment)
66
	{
67
		lwline_free(root->segment);
68
	}
69 70
	lwfree(root);
}
71

72 73 74
/**
* Free the cache object and all the sub-objects properly.
*/
75
static void
76 77 78 79 80 81
RTreeCacheClear(RTREE_POLY_CACHE* cache)
{
	int g, r, i;
	POSTGIS_DEBUGF(2, "RTreeCacheClear called for %p", cache);
	i = 0;
	for (g = 0; g < cache->polyCount; g++)
82
	{
83
		for (r = 0; r < cache->ringCounts[g]; r++)
84
		{
85
			RTreeFree(cache->ringIndices[i]);
86 87 88
			i++;
		}
	}
89 90 91 92 93 94
	lwfree(cache->ringIndices);
	lwfree(cache->ringCounts);
	cache->ringIndices = 0;
	cache->ringCounts = 0;
	cache->polyCount = 0;
}
95 96


97
/**
98
 * Returns 1 if min < value <= max, 0 otherwise.
99
*/
100
static uint32
101 102 103 104 105 106 107 108
IntervalIsContained(RTREE_INTERVAL* interval, double value)
{
	return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
}

/**
* Creates an interval with the total extents of the two given intervals.
*/
109
static RTREE_INTERVAL*
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
RTreeMergeIntervals(RTREE_INTERVAL *inter1, RTREE_INTERVAL *inter2)
{
	RTREE_INTERVAL *interval;

	POSTGIS_DEBUGF(2, "RTreeMergeIntervals called with %p, %p", inter1, inter2);

	interval = lwalloc(sizeof(RTREE_INTERVAL));
	interval->max = FP_MAX(inter1->max, inter2->max);
	interval->min = FP_MIN(inter1->min, inter2->min);

	POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);

	return interval;
}

/**
* Creates an interval given the min and max values, in arbitrary order.
*/
static RTREE_INTERVAL*
RTreeCreateInterval(double value1, double value2)
{
	RTREE_INTERVAL *interval;

	POSTGIS_DEBUGF(2, "RTreeCreateInterval called with %8.3f, %8.3f", value1, value2);

	interval = lwalloc(sizeof(RTREE_INTERVAL));
	interval->max = FP_MAX(value1, value2);
	interval->min = FP_MIN(value1, value2);

	POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);

	return interval;
142 143 144
}

/**
145 146
* Creates an interior node given the children.
*/
147
static RTREE_NODE*
148
RTreeCreateInteriorNode(RTREE_NODE* left, RTREE_NODE* right)
149 150 151
{
	RTREE_NODE *parent;

152
	POSTGIS_DEBUGF(2, "RTreeCreateInteriorNode called for children %p, %p", left, right);
153 154 155 156

	parent = lwalloc(sizeof(RTREE_NODE));
	parent->leftNode = left;
	parent->rightNode = right;
157
	parent->interval = RTreeMergeIntervals(left->interval, right->interval);
158 159
	parent->segment = NULL;

160
	POSTGIS_DEBUGF(3, "RTreeCreateInteriorNode returning %p", parent);
161 162 163 164 165

	return parent;
}

/**
166 167
* Creates a leaf node given the pointer to the start point of the segment.
*/
168
static RTREE_NODE*
169
RTreeCreateLeafNode(POINTARRAY* pa, int startPoint)
170 171 172 173 174 175 176 177
{
	RTREE_NODE *parent;
	LWLINE *line;
	double value1;
	double value2;
	POINT4D tmp;
	POINTARRAY *npa;

178
	POSTGIS_DEBUGF(2, "RTreeCreateLeafNode called for point %d of %p", startPoint, pa);
179 180 181

	if (pa->npoints < startPoint + 2)
	{
182
		lwpgerror("RTreeCreateLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
183 184 185 186 187 188 189
	}

	/*
	 * 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.
	 */
190
	npa = ptarray_construct_empty(0,0,2);
191 192 193

	getPoint4d_p(pa, startPoint, &tmp);
	value1 = tmp.y;
194 195 196
	ptarray_append_point(npa,&tmp,LW_TRUE);
	
	getPoint4d_p(pa, startPoint+1, &tmp);
197
	value2 = tmp.y;
198
	ptarray_append_point(npa,&tmp,LW_TRUE);
199

200 201
	line = lwline_construct(SRID_UNKNOWN, NULL, npa);
	
202
	parent = lwalloc(sizeof(RTREE_NODE));
203
	parent->interval = RTreeCreateInterval(value1, value2);
204 205 206 207
	parent->segment = line;
	parent->leftNode = NULL;
	parent->rightNode = NULL;

208
	POSTGIS_DEBUGF(3, "RTreeCreateLeafNode returning %p", parent);
209 210 211 212 213

	return parent;
}

/**
214 215 216
* Creates an rtree given a pointer to the point array.
* Must copy the point array.
*/
217
static RTREE_NODE*
218
RTreeCreate(POINTARRAY* pointArray)
219
{
220 221 222 223
	RTREE_NODE* root;
	RTREE_NODE** nodes = lwalloc(pointArray->npoints * sizeof(RTREE_NODE*));
	int i, nodeCount;
	int childNodes, parentNodes;
224

225
	POSTGIS_DEBUGF(2, "RTreeCreate called with pointarray %p", pointArray);
226

227
	nodeCount = pointArray->npoints - 1;
228

229
	POSTGIS_DEBUGF(3, "Total leaf nodes: %d", nodeCount);
230

231 232 233 234
	/*
	 * Create a leaf node for every line segment.
	 */
	for (i = 0; i < nodeCount; i++)
235
	{
236
		nodes[i] = RTreeCreateLeafNode(pointArray, i);
237 238
	}

239 240 241 242 243 244 245 246
	/*
	 * 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)
247
	{
248
		POSTGIS_DEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes);
249

250 251
		i = 0;
		while (i < parentNodes)
252
		{
253 254
			nodes[i] = RTreeCreateInteriorNode(nodes[i*2], nodes[i*2+1]);
			i++;
255
		}
256 257 258 259
		/*
		 * Check for an odd numbered final node.
		 */
		if (parentNodes * 2 < childNodes)
260
		{
261
			POSTGIS_DEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i);
262

263 264
			nodes[i] = nodes[childNodes - 1];
			parentNodes++;
265
		}
266 267
		childNodes = parentNodes;
		parentNodes = parentNodes / 2;
268 269
	}

270 271 272 273 274
	root = nodes[0];
	lwfree(nodes);
	POSTGIS_DEBUGF(3, "RTreeCreate returning %p", root);

	return root;
275 276
}

277

278 279
/**
* Merges two multilinestrings into a single multilinestring.
280
*/
281
static LWMLINE*
282
RTreeMergeMultiLines(LWMLINE *line1, LWMLINE *line2)
283 284 285 286 287
{
	LWGEOM **geoms;
	LWCOLLECTION *col;
	int i, j, ngeoms;

288
	POSTGIS_DEBUGF(2, "RTreeMergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, line1->type, line2, line2->ngeoms, line2->type);
289 290 291 292 293 294 295 296 297 298 299 300 301

	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]);
	}
302
	col = lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, ngeoms, geoms);
303

304
	POSTGIS_DEBUGF(3, "RTreeMergeMultiLines returning %p, %d, %d", col, col->ngeoms, col->type);
305 306 307 308 309

	return (LWMLINE *)col;
}


310 311 312 313 314
/**
* Callback function sent into the GetGeomCache generic caching system. Given a
* LWGEOM* this function builds and stores an RTREE_POLY_CACHE into the provided
* GeomCache object.
*/
315
static int
316
RTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
317
{
318
	int i, p, r;
319 320 321
	LWMPOLY *mpoly;
	LWPOLY *poly;
	int nrings;
322 323 324 325 326
	RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
	RTREE_POLY_CACHE* currentCache;
	
	if ( ! cache )
		return LW_FAILURE;
327

328 329
	if ( rtree_cache->index )
	{
330
		lwpgerror("RTreeBuilder asked to build index where one already exists.");
331 332 333
		return LW_FAILURE;
	}
	
334
	if (lwgeom->type == MULTIPOLYGONTYPE)
335
	{
336
		POSTGIS_DEBUG(2, "RTreeBuilder MULTIPOLYGON");
337 338 339 340 341
		mpoly = (LWMPOLY *)lwgeom;
		nrings = 0;
		/*
		** Count the total number of rings.
		*/
342
		currentCache = RTreeCacheCreate();
343 344
		currentCache->polyCount = mpoly->ngeoms;
		currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms);
345 346
		for ( i = 0; i < mpoly->ngeoms; i++ )
		{
347
			currentCache->ringCounts[i] = mpoly->geoms[i]->nrings;
348 349 350 351
			nrings += mpoly->geoms[i]->nrings;
		}
		currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings);
		/*
352 353
		** Load the array in geometry order, each outer ring followed by the inner rings
                ** associated with that outer ring
354
		*/
355 356
                i = 0;
		for ( p = 0; p < mpoly->ngeoms; p++ )
357
		{
358
			for ( r = 0; r < mpoly->geoms[p]->nrings; r++ )
359
			{
360
				currentCache->ringIndices[i] = RTreeCreate(mpoly->geoms[p]->rings[r]);
361 362 363
				i++;
			}
		}
364
		rtree_cache->index = currentCache;
365
	}
366
	else if ( lwgeom->type == POLYGONTYPE )
367
	{
368
		POSTGIS_DEBUG(2, "RTreeBuilder POLYGON");
369
		poly = (LWPOLY *)lwgeom;
370
		currentCache = RTreeCacheCreate();
371
		currentCache->polyCount = 1;
372 373
		currentCache->ringCounts = lwalloc(sizeof(int));
		currentCache->ringCounts[0] = poly->nrings;
374 375 376 377 378 379
		/*
		** Just load the rings on in order
		*/
		currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
		for ( i = 0; i < poly->nrings; i++ )
		{
380
			currentCache->ringIndices[i] = RTreeCreate(poly->rings[i]);
381
		}
382
		rtree_cache->index = currentCache;
383 384 385 386
	}
	else
	{
		/* Uh oh, shouldn't be here. */
387
		lwpgerror("RTreeBuilder got asked to build index on non-polygon");
388
		return LW_FAILURE;
389
	}
390 391
	return LW_SUCCESS;	
}
392

393
/**
394
* Callback function sent into the GetGeomCache generic caching system. On a
395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
* cache miss, this function clears the cached index object.
*/
static int
RTreeFreer(GeomCache* cache)
{
	RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
	
	if ( ! cache )
		return LW_FAILURE;
	
	if ( rtree_cache->index )
	{
		RTreeCacheClear(rtree_cache->index);
		lwfree(rtree_cache->index);
		rtree_cache->index = 0;
		rtree_cache->argnum = 0;
	}
	return LW_SUCCESS;
}

static GeomCache*
RTreeAllocator(void)
{
	RTreeGeomCache* cache = palloc(sizeof(RTreeGeomCache));
	memset(cache, 0, sizeof(RTreeGeomCache));
	return (GeomCache*)cache;
}

static GeomCacheMethods RTreeCacheMethods =
{
	RTREE_CACHE_ENTRY,
	RTreeBuilder,
	RTreeFreer,
	RTreeAllocator
};

RTREE_POLY_CACHE*
GetRtreeCache(FunctionCallInfoData* fcinfo, GSERIALIZED* g1)
{
	RTreeGeomCache* cache = (RTreeGeomCache*)GetGeomCache(fcinfo, &RTreeCacheMethods, g1, NULL);
	RTREE_POLY_CACHE* index = NULL;

	if ( cache )
		index = cache->index;

	return index;
441 442
}

443

444
/**
445 446 447 448 449 450
* Retrieves a collection of line segments given the root and crossing value.
* The collection is a multilinestring consisting of two point lines
* representing the segments of the ring that may be crossed by the
* horizontal projection line at the given y value.
*/
LWMLINE *RTreeFindLineSegments(RTREE_NODE *root, double value)
451
{
452 453
	LWMLINE *tmp, *result;
	LWGEOM **lwgeoms;
454

455
	POSTGIS_DEBUGF(2, "RTreeFindLineSegments called for tree %p and value %8.3f", root, value);
456

457
	result = NULL;
458

459
	if (!IntervalIsContained(root->interval, value))
460
	{
461 462 463
		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: not contained.", root);

		return NULL;
464
	}
465 466 467

	/* If there is a segment defined for this node, include it. */
	if (root->segment)
468
	{
469 470 471 472
		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: adding segment %p %d.", root, root->segment, root->segment->type);

		lwgeoms = lwalloc(sizeof(LWGEOM *));
		lwgeoms[0] = (LWGEOM *)root->segment;
473

474
		POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", root->segment, root->segment->type, FLAGS_GET_Z(root->segment->flags));
475

476
		result = (LWMLINE *)lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, 1, lwgeoms);
477
	}
478 479 480

	/* If there is a left child node, recursively include its results. */
	if (root->leftNode)
481
	{
482 483 484 485 486 487 488 489 490 491 492 493
		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing left.", root);

		tmp = RTreeFindLineSegments(root->leftNode, value);
		if (tmp)
		{
			POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));

			if (result)
				result = RTreeMergeMultiLines(result, tmp);
			else
				result = tmp;
		}
494 495
	}

496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
	/* Same for any right child. */
	if (root->rightNode)
	{
		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing right.", root);

		tmp = RTreeFindLineSegments(root->rightNode, value);
		if (tmp)
		{
			POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));

			if (result)
				result = RTreeMergeMultiLines(result, tmp);
			else
				result = tmp;
		}
	}
512

513
	return result;
514 515
}

516 517