PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwgeom_rtree.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  * Copyright (C) 2001-2005 Refractions Research Inc.
6  *
7  * This is free software; you can redistribute and/or modify it under
8  * the terms of the GNU General Public Licence. See the COPYING file.
9  *
10  **********************************************************************/
11 
12 #include <assert.h>
13 
14 #include "../postgis_config.h"
15 #include "lwgeom_pg.h"
16 #include "liblwgeom.h"
17 #include "liblwgeom_internal.h" /* For FP comparators. */
18 #include "lwgeom_cache.h"
19 #include "lwgeom_rtree.h"
20 
21 
22 /* Prototypes */
23 static void RTreeFree(RTREE_NODE* root);
24 
28 static RTREE_POLY_CACHE*
30 {
32  result = lwalloc(sizeof(RTREE_POLY_CACHE));
33  memset(result, 0, sizeof(RTREE_POLY_CACHE));
34  return result;
35 }
36 
41 static void
43 {
44  POSTGIS_DEBUGF(2, "RTreeFree called for %p", root);
45 
46  if (root->leftNode)
47  RTreeFree(root->leftNode);
48  if (root->rightNode)
49  RTreeFree(root->rightNode);
50  lwfree(root->interval);
51  if (root->segment)
52  {
53  lwline_free(root->segment);
54  }
55  lwfree(root);
56 }
57 
61 static void
63 {
64  int g, r, i;
65  POSTGIS_DEBUGF(2, "RTreeCacheClear called for %p", cache);
66  i = 0;
67  for (g = 0; g < cache->polyCount; g++)
68  {
69  for (r = 0; r < cache->ringCounts[g]; r++)
70  {
71  RTreeFree(cache->ringIndices[i]);
72  i++;
73  }
74  }
75  lwfree(cache->ringIndices);
76  lwfree(cache->ringCounts);
77  cache->ringIndices = 0;
78  cache->ringCounts = 0;
79  cache->polyCount = 0;
80 }
81 
82 
86 static uint32
88 {
89  return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
90 }
91 
95 static RTREE_INTERVAL*
97 {
98  RTREE_INTERVAL *interval;
99 
100  POSTGIS_DEBUGF(2, "RTreeMergeIntervals called with %p, %p", inter1, inter2);
101 
102  interval = lwalloc(sizeof(RTREE_INTERVAL));
103  interval->max = FP_MAX(inter1->max, inter2->max);
104  interval->min = FP_MIN(inter1->min, inter2->min);
105 
106  POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);
107 
108  return interval;
109 }
110 
114 static RTREE_INTERVAL*
115 RTreeCreateInterval(double value1, double value2)
116 {
117  RTREE_INTERVAL *interval;
118 
119  POSTGIS_DEBUGF(2, "RTreeCreateInterval called with %8.3f, %8.3f", value1, value2);
120 
121  interval = lwalloc(sizeof(RTREE_INTERVAL));
122  interval->max = FP_MAX(value1, value2);
123  interval->min = FP_MIN(value1, value2);
124 
125  POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);
126 
127  return interval;
128 }
129 
133 static RTREE_NODE*
135 {
136  RTREE_NODE *parent;
137 
138  POSTGIS_DEBUGF(2, "RTreeCreateInteriorNode called for children %p, %p", left, right);
139 
140  parent = lwalloc(sizeof(RTREE_NODE));
141  parent->leftNode = left;
142  parent->rightNode = right;
143  parent->interval = RTreeMergeIntervals(left->interval, right->interval);
144  parent->segment = NULL;
145 
146  POSTGIS_DEBUGF(3, "RTreeCreateInteriorNode returning %p", parent);
147 
148  return parent;
149 }
150 
154 static RTREE_NODE*
155 RTreeCreateLeafNode(POINTARRAY* pa, int startPoint)
156 {
157  RTREE_NODE *parent;
158  LWLINE *line;
159  double value1;
160  double value2;
161  POINT4D tmp;
162  POINTARRAY *npa;
163 
164  POSTGIS_DEBUGF(2, "RTreeCreateLeafNode called for point %d of %p", startPoint, pa);
165 
166  if (pa->npoints < startPoint + 2)
167  {
168  lwerror("RTreeCreateLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
169  }
170 
171  /*
172  * The given point array will be part of a geometry that will be freed
173  * independently of the index. Since we may want to cache the index,
174  * we must create independent arrays.
175  */
176  npa = ptarray_construct_empty(0,0,2);
177 
178  getPoint4d_p(pa, startPoint, &tmp);
179  value1 = tmp.y;
180  ptarray_append_point(npa,&tmp,LW_TRUE);
181 
182  getPoint4d_p(pa, startPoint+1, &tmp);
183  value2 = tmp.y;
184  ptarray_append_point(npa,&tmp,LW_TRUE);
185 
186  line = lwline_construct(SRID_UNKNOWN, NULL, npa);
187 
188  parent = lwalloc(sizeof(RTREE_NODE));
189  parent->interval = RTreeCreateInterval(value1, value2);
190  parent->segment = line;
191  parent->leftNode = NULL;
192  parent->rightNode = NULL;
193 
194  POSTGIS_DEBUGF(3, "RTreeCreateLeafNode returning %p", parent);
195 
196  return parent;
197 }
198 
203 static RTREE_NODE*
205 {
206  RTREE_NODE* root;
207  RTREE_NODE** nodes = lwalloc(pointArray->npoints * sizeof(RTREE_NODE*));
208  int i, nodeCount;
209  int childNodes, parentNodes;
210 
211  POSTGIS_DEBUGF(2, "RTreeCreate called with pointarray %p", pointArray);
212 
213  nodeCount = pointArray->npoints - 1;
214 
215  POSTGIS_DEBUGF(3, "Total leaf nodes: %d", nodeCount);
216 
217  /*
218  * Create a leaf node for every line segment.
219  */
220  for (i = 0; i < nodeCount; i++)
221  {
222  nodes[i] = RTreeCreateLeafNode(pointArray, i);
223  }
224 
225  /*
226  * Next we group nodes by pairs. If there's an odd number of nodes,
227  * we bring the last node up a level as is. Continue until we have
228  * a single top node.
229  */
230  childNodes = nodeCount;
231  parentNodes = nodeCount / 2;
232  while (parentNodes > 0)
233  {
234  POSTGIS_DEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes);
235 
236  i = 0;
237  while (i < parentNodes)
238  {
239  nodes[i] = RTreeCreateInteriorNode(nodes[i*2], nodes[i*2+1]);
240  i++;
241  }
242  /*
243  * Check for an odd numbered final node.
244  */
245  if (parentNodes * 2 < childNodes)
246  {
247  POSTGIS_DEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i);
248 
249  nodes[i] = nodes[childNodes - 1];
250  parentNodes++;
251  }
252  childNodes = parentNodes;
253  parentNodes = parentNodes / 2;
254  }
255 
256  root = nodes[0];
257  lwfree(nodes);
258  POSTGIS_DEBUGF(3, "RTreeCreate returning %p", root);
259 
260  return root;
261 }
262 
263 
267 static LWMLINE*
269 {
270  LWGEOM **geoms;
271  LWCOLLECTION *col;
272  int i, j, ngeoms;
273 
274  POSTGIS_DEBUGF(2, "RTreeMergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, line1->type, line2, line2->ngeoms, line2->type);
275 
276  ngeoms = line1->ngeoms + line2->ngeoms;
277  geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);
278 
279  j = 0;
280  for (i = 0; i < line1->ngeoms; i++, j++)
281  {
282  geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);
283  }
284  for (i = 0; i < line2->ngeoms; i++, j++)
285  {
286  geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);
287  }
288  col = lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, ngeoms, geoms);
289 
290  POSTGIS_DEBUGF(3, "RTreeMergeMultiLines returning %p, %d, %d", col, col->ngeoms, col->type);
291 
292  return (LWMLINE *)col;
293 }
294 
295 
301 static int
302 RTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
303 {
304  int i, p, r;
305  LWMPOLY *mpoly;
306  LWPOLY *poly;
307  int nrings;
308  RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
309  RTREE_POLY_CACHE* currentCache;
310 
311  if ( ! cache )
312  return LW_FAILURE;
313 
314  if ( rtree_cache->index )
315  {
316  lwerror("RTreeBuilder asked to build index where one already exists.");
317  return LW_FAILURE;
318  }
319 
320  if (lwgeom->type == MULTIPOLYGONTYPE)
321  {
322  POSTGIS_DEBUG(2, "RTreeBuilder MULTIPOLYGON");
323  mpoly = (LWMPOLY *)lwgeom;
324  nrings = 0;
325  /*
326  ** Count the total number of rings.
327  */
328  currentCache = RTreeCacheCreate();
329  currentCache->polyCount = mpoly->ngeoms;
330  currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms);
331  for ( i = 0; i < mpoly->ngeoms; i++ )
332  {
333  currentCache->ringCounts[i] = mpoly->geoms[i]->nrings;
334  nrings += mpoly->geoms[i]->nrings;
335  }
336  currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings);
337  /*
338  ** Load the array in geometry order, each outer ring followed by the inner rings
339  ** associated with that outer ring
340  */
341  i = 0;
342  for ( p = 0; p < mpoly->ngeoms; p++ )
343  {
344  for ( r = 0; r < mpoly->geoms[p]->nrings; r++ )
345  {
346  currentCache->ringIndices[i] = RTreeCreate(mpoly->geoms[p]->rings[r]);
347  i++;
348  }
349  }
350  rtree_cache->index = currentCache;
351  }
352  else if ( lwgeom->type == POLYGONTYPE )
353  {
354  POSTGIS_DEBUG(2, "RTreeBuilder POLYGON");
355  poly = (LWPOLY *)lwgeom;
356  currentCache = RTreeCacheCreate();
357  currentCache->polyCount = 1;
358  currentCache->ringCounts = lwalloc(sizeof(int));
359  currentCache->ringCounts[0] = poly->nrings;
360  /*
361  ** Just load the rings on in order
362  */
363  currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
364  for ( i = 0; i < poly->nrings; i++ )
365  {
366  currentCache->ringIndices[i] = RTreeCreate(poly->rings[i]);
367  }
368  rtree_cache->index = currentCache;
369  }
370  else
371  {
372  /* Uh oh, shouldn't be here. */
373  lwerror("RTreeBuilder got asked to build index on non-polygon");
374  return LW_FAILURE;
375  }
376  return LW_SUCCESS;
377 }
378 
383 static int
384 RTreeFreer(GeomCache* cache)
385 {
386  RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
387 
388  if ( ! cache )
389  return LW_FAILURE;
390 
391  if ( rtree_cache->index )
392  {
393  RTreeCacheClear(rtree_cache->index);
394  lwfree(rtree_cache->index);
395  rtree_cache->index = 0;
396  rtree_cache->argnum = 0;
397  }
398  return LW_SUCCESS;
399 }
400 
401 static GeomCache*
403 {
404  RTreeGeomCache* cache = palloc(sizeof(RTreeGeomCache));
405  memset(cache, 0, sizeof(RTreeGeomCache));
406  return (GeomCache*)cache;
407 }
408 
409 static GeomCacheMethods RTreeCacheMethods =
410 {
411  RTREE_CACHE_ENTRY,
412  RTreeBuilder,
413  RTreeFreer,
415 };
416 
418 GetRtreeCache(FunctionCallInfoData* fcinfo, GSERIALIZED* g1)
419 {
420  RTreeGeomCache* cache = (RTreeGeomCache*)GetGeomCache(fcinfo, &RTreeCacheMethods, g1, NULL);
421  RTREE_POLY_CACHE* index = NULL;
422 
423  if ( cache )
424  index = cache->index;
425 
426  return index;
427 }
428 
429 
437 {
438  LWMLINE *tmp, *result;
439  LWGEOM **lwgeoms;
440 
441  POSTGIS_DEBUGF(2, "RTreeFindLineSegments called for tree %p and value %8.3f", root, value);
442 
443  result = NULL;
444 
445  if (!IntervalIsContained(root->interval, value))
446  {
447  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: not contained.", root);
448 
449  return NULL;
450  }
451 
452  /* If there is a segment defined for this node, include it. */
453  if (root->segment)
454  {
455  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: adding segment %p %d.", root, root->segment, root->segment->type);
456 
457  lwgeoms = lwalloc(sizeof(LWGEOM *));
458  lwgeoms[0] = (LWGEOM *)root->segment;
459 
460  POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", root->segment, root->segment->type, FLAGS_GET_Z(root->segment->flags));
461 
462  result = (LWMLINE *)lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, 1, lwgeoms);
463  }
464 
465  /* If there is a left child node, recursively include its results. */
466  if (root->leftNode)
467  {
468  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing left.", root);
469 
470  tmp = RTreeFindLineSegments(root->leftNode, value);
471  if (tmp)
472  {
473  POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));
474 
475  if (result)
476  result = RTreeMergeMultiLines(result, tmp);
477  else
478  result = tmp;
479  }
480  }
481 
482  /* Same for any right child. */
483  if (root->rightNode)
484  {
485  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing right.", root);
486 
487  tmp = RTreeFindLineSegments(root->rightNode, value);
488  if (tmp)
489  {
490  POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));
491 
492  if (result)
493  result = RTreeMergeMultiLines(result, tmp);
494  else
495  result = tmp;
496  }
497  }
498 
499  return result;
500 }
501 
502 
503 
static int RTreeFreer(GeomCache *cache)
Callback function sent into the GetGeomCache generic caching system.
Definition: lwgeom_rtree.c:384
uint8_t type
Definition: liblwgeom.h:374
uint8_t type
Definition: liblwgeom.h:433
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
char * r
Definition: cu_in_wkt.c:25
static LWMLINE * RTreeMergeMultiLines(LWMLINE *line1, LWMLINE *line2)
Merges two multilinestrings into a single multilinestring.
Definition: lwgeom_rtree.c:268
void lwfree(void *mem)
Definition: lwutil.c:190
int npoints
Definition: liblwgeom.h:327
uint8_t type
Definition: liblwgeom.h:459
static GeomCache * RTreeAllocator(void)
Definition: lwgeom_rtree.c:402
The following struct and methods are used for a 1D RTree implementation, described at: http://lin-ear...
Definition: lwgeom_rtree.h:21
#define POLYGONTYPE
Definition: liblwgeom.h:62
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:57
The tree structure used for fast P-i-P tests by point_in_multipolygon_rtree()
Definition: lwgeom_rtree.h:33
#define LW_SUCCESS
Definition: liblwgeom.h:55
void lwline_free(LWLINE *line)
Definition: lwline.c:63
struct rtree_node * leftNode
Definition: lwgeom_rtree.h:24
LWMLINE * RTreeFindLineSegments(RTREE_NODE *root, double value)
Retrieves a collection of line segments given the root and crossing value.
Definition: lwgeom_rtree.c:436
static void RTreeFree(RTREE_NODE *root)
Recursively frees the child nodes, the interval and the line before freeing the root node...
Definition: lwgeom_rtree.c:42
#define FP_MIN(A, B)
Representation for the y-axis interval spanned by an edge.
Definition: lwgeom_rtree.h:9
char ** result
Definition: liblwgeom.h:218
int ngeoms
Definition: liblwgeom.h:437
#define FP_CONTAINS_INCL(A, X, B)
static uint32 IntervalIsContained(RTREE_INTERVAL *interval, double value)
Returns 1 if min < value <= max, 0 otherwise.
Definition: lwgeom_rtree.c:87
#define LW_FAILURE
Definition: liblwgeom.h:54
static RTREE_NODE * RTreeCreateInteriorNode(RTREE_NODE *left, RTREE_NODE *right)
Creates an interior node given the children.
Definition: lwgeom_rtree.c:134
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
RTREE_INTERVAL * interval
Definition: lwgeom_rtree.h:23
uint8_t flags
Definition: liblwgeom.h:434
unsigned int uint32
Definition: shpopen.c:274
static RTREE_INTERVAL * RTreeMergeIntervals(RTREE_INTERVAL *inter1, RTREE_INTERVAL *inter2)
Creates an interval with the total extents of the two given intervals.
Definition: lwgeom_rtree.c:96
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_TRUE, then a duplicate point will not be added.
Definition: ptarray.c:141
static RTREE_NODE * RTreeCreate(POINTARRAY *pointArray)
Creates an rtree given a pointer to the point array.
Definition: lwgeom_rtree.c:204
LWPOLY ** geoms
Definition: liblwgeom.h:452
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
static RTREE_POLY_CACHE * RTreeCacheCreate()
Allocate a fresh clean RTREE_POLY_CACHE.
Definition: lwgeom_rtree.c:29
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:154
static RTREE_NODE * RTreeCreateLeafNode(POINTARRAY *pa, int startPoint)
Creates a leaf node given the pointer to the start point of the segment.
Definition: lwgeom_rtree.c:155
POINTARRAY ** rings
Definition: liblwgeom.h:413
int nrings
Definition: liblwgeom.h:411
static int RTreeBuilder(const LWGEOM *lwgeom, GeomCache *cache)
Callback function sent into the GetGeomCache generic caching system.
Definition: lwgeom_rtree.c:302
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
static GeomCacheMethods RTreeCacheMethods
Definition: lwgeom_rtree.c:409
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:351
int ngeoms
Definition: liblwgeom.h:450
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
LWLINE * segment
Definition: lwgeom_rtree.h:26
LWLINE ** geoms
Definition: liblwgeom.h:439
static void RTreeCacheClear(RTREE_POLY_CACHE *cache)
Free the cache object and all the sub-objects properly.
Definition: lwgeom_rtree.c:62
struct rtree_node * rightNode
Definition: lwgeom_rtree.h:25
static RTREE_INTERVAL * RTreeCreateInterval(double value1, double value2)
Creates an interval given the min and max values, in arbitrary order.
Definition: lwgeom_rtree.c:115
uint8_t type
Definition: liblwgeom.h:352
void * lwalloc(size_t size)
Definition: lwutil.c:175
RTREE_NODE ** ringIndices
Definition: lwgeom_rtree.h:35
RTREE_POLY_CACHE * GetRtreeCache(FunctionCallInfoData *fcinfo, GSERIALIZED *g1)
Checks for a cache hit against the provided geometry and returns a pre-built index structure (RTREE_P...
Definition: lwgeom_rtree.c:418
double y
Definition: liblwgeom.h:308
#define MULTILINETYPE
Definition: liblwgeom.h:64
uint8_t flags
Definition: liblwgeom.h:375
RTREE_POLY_CACHE * index
Definition: lwgeom_rtree.h:50
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:217
This library is the generic geometry handling section of PostGIS.
#define FP_MAX(A, B)