PostGIS  2.5.0dev-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  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2001-2005 Refractions Research Inc.
22  *
23  **********************************************************************/
24 
25 
26 #include <assert.h>
27 
28 #include "../postgis_config.h"
29 #include "lwgeom_pg.h"
30 #include "liblwgeom.h"
31 #include "liblwgeom_internal.h" /* For FP comparators. */
32 #include "lwgeom_cache.h"
33 #include "lwgeom_rtree.h"
34 
35 
36 /* Prototypes */
37 static void RTreeFree(RTREE_NODE* root);
38 
42 static RTREE_POLY_CACHE*
44 {
45  RTREE_POLY_CACHE* result;
46  result = lwalloc(sizeof(RTREE_POLY_CACHE));
47  memset(result, 0, sizeof(RTREE_POLY_CACHE));
48  return result;
49 }
50 
55 static void
57 {
58  POSTGIS_DEBUGF(2, "RTreeFree called for %p", root);
59 
60  if (root->leftNode)
61  RTreeFree(root->leftNode);
62  if (root->rightNode)
63  RTreeFree(root->rightNode);
64  lwfree(root->interval);
65  if (root->segment)
66  {
67  lwline_free(root->segment);
68  }
69  lwfree(root);
70 }
71 
75 static void
77 {
78  int g, r, i;
79  POSTGIS_DEBUGF(2, "RTreeCacheClear called for %p", cache);
80  i = 0;
81  for (g = 0; g < cache->polyCount; g++)
82  {
83  for (r = 0; r < cache->ringCounts[g]; r++)
84  {
85  RTreeFree(cache->ringIndices[i]);
86  i++;
87  }
88  }
89  lwfree(cache->ringIndices);
90  lwfree(cache->ringCounts);
91  cache->ringIndices = 0;
92  cache->ringCounts = 0;
93  cache->polyCount = 0;
94 }
95 
96 
100 static uint32
102 {
103  return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
104 }
105 
109 static RTREE_INTERVAL*
111 {
112  RTREE_INTERVAL *interval;
113 
114  POSTGIS_DEBUGF(2, "RTreeMergeIntervals called with %p, %p", inter1, inter2);
115 
116  interval = lwalloc(sizeof(RTREE_INTERVAL));
117  interval->max = FP_MAX(inter1->max, inter2->max);
118  interval->min = FP_MIN(inter1->min, inter2->min);
119 
120  POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);
121 
122  return interval;
123 }
124 
128 static RTREE_INTERVAL*
129 RTreeCreateInterval(double value1, double value2)
130 {
131  RTREE_INTERVAL *interval;
132 
133  POSTGIS_DEBUGF(2, "RTreeCreateInterval called with %8.3f, %8.3f", value1, value2);
134 
135  interval = lwalloc(sizeof(RTREE_INTERVAL));
136  interval->max = FP_MAX(value1, value2);
137  interval->min = FP_MIN(value1, value2);
138 
139  POSTGIS_DEBUGF(3, "interval min = %8.3f, max = %8.3f", interval->min, interval->max);
140 
141  return interval;
142 }
143 
147 static RTREE_NODE*
149 {
150  RTREE_NODE *parent;
151 
152  POSTGIS_DEBUGF(2, "RTreeCreateInteriorNode called for children %p, %p", left, right);
153 
154  parent = lwalloc(sizeof(RTREE_NODE));
155  parent->leftNode = left;
156  parent->rightNode = right;
157  parent->interval = RTreeMergeIntervals(left->interval, right->interval);
158  parent->segment = NULL;
159 
160  POSTGIS_DEBUGF(3, "RTreeCreateInteriorNode returning %p", parent);
161 
162  return parent;
163 }
164 
168 static RTREE_NODE*
170 {
171  RTREE_NODE *parent;
172  LWLINE *line;
173  double value1;
174  double value2;
175  POINT4D tmp;
176  POINTARRAY *npa;
177 
178  POSTGIS_DEBUGF(2, "RTreeCreateLeafNode called for point %d of %p", startPoint, pa);
179 
180  if (pa->npoints < startPoint + 2)
181  {
182  lwpgerror("RTreeCreateLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
183  }
184 
185  /*
186  * The given point array will be part of a geometry that will be freed
187  * independently of the index. Since we may want to cache the index,
188  * we must create independent arrays.
189  */
190  npa = ptarray_construct_empty(0,0,2);
191 
192  getPoint4d_p(pa, startPoint, &tmp);
193  value1 = tmp.y;
194  ptarray_append_point(npa,&tmp,LW_TRUE);
195 
196  getPoint4d_p(pa, startPoint+1, &tmp);
197  value2 = tmp.y;
198  ptarray_append_point(npa,&tmp,LW_TRUE);
199 
200  line = lwline_construct(SRID_UNKNOWN, NULL, npa);
201 
202  parent = lwalloc(sizeof(RTREE_NODE));
203  parent->interval = RTreeCreateInterval(value1, value2);
204  parent->segment = line;
205  parent->leftNode = NULL;
206  parent->rightNode = NULL;
207 
208  POSTGIS_DEBUGF(3, "RTreeCreateLeafNode returning %p", parent);
209 
210  return parent;
211 }
212 
217 static RTREE_NODE*
219 {
220  RTREE_NODE* root;
221  RTREE_NODE** nodes = lwalloc(pointArray->npoints * sizeof(RTREE_NODE*));
222  uint32_t i, nodeCount;
223  uint32_t 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  * Create a leaf node for every line segment.
233  */
234  for (i = 0; i < nodeCount; i++)
235  {
236  nodes[i] = RTreeCreateLeafNode(pointArray, i);
237  }
238 
239  /*
240  * Next we group nodes by pairs. If there's an odd number of nodes,
241  * we bring the last node up a level as is. Continue until we have
242  * a single top node.
243  */
244  childNodes = nodeCount;
245  parentNodes = nodeCount / 2;
246  while (parentNodes > 0)
247  {
248  POSTGIS_DEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes);
249 
250  i = 0;
251  while (i < parentNodes)
252  {
253  nodes[i] = RTreeCreateInteriorNode(nodes[i*2], nodes[i*2+1]);
254  i++;
255  }
256  /*
257  * Check for an odd numbered final node.
258  */
259  if (parentNodes * 2 < childNodes)
260  {
261  POSTGIS_DEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i);
262 
263  nodes[i] = nodes[childNodes - 1];
264  parentNodes++;
265  }
266  childNodes = parentNodes;
267  parentNodes = parentNodes / 2;
268  }
269 
270  root = nodes[0];
271  lwfree(nodes);
272  POSTGIS_DEBUGF(3, "RTreeCreate returning %p", root);
273 
274  return root;
275 }
276 
277 
281 static LWMLINE*
283 {
284  LWGEOM **geoms;
285  LWCOLLECTION *col;
286  uint32_t i, j, ngeoms;
287 
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  ngeoms = line1->ngeoms + line2->ngeoms;
291  geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);
292 
293  j = 0;
294  for (i = 0; i < line1->ngeoms; i++, j++)
295  {
296  geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);
297  }
298  for (i = 0; i < line2->ngeoms; i++, j++)
299  {
300  geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);
301  }
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  return (LWMLINE *)col;
307 }
308 
309 
315 static int
316 RTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
317 {
318  uint32_t i, p, r;
319  LWMPOLY *mpoly;
320  LWPOLY *poly;
321  int nrings;
322  RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
323  RTREE_POLY_CACHE* currentCache;
324 
325  if ( ! cache )
326  return LW_FAILURE;
327 
328  if ( rtree_cache->index )
329  {
330  lwpgerror("RTreeBuilder asked to build index where one already exists.");
331  return LW_FAILURE;
332  }
333 
334  if (lwgeom->type == MULTIPOLYGONTYPE)
335  {
336  POSTGIS_DEBUG(2, "RTreeBuilder MULTIPOLYGON");
337  mpoly = (LWMPOLY *)lwgeom;
338  nrings = 0;
339  /*
340  ** Count the total number of rings.
341  */
342  currentCache = RTreeCacheCreate();
343  currentCache->polyCount = mpoly->ngeoms;
344  currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms);
345  for ( i = 0; i < mpoly->ngeoms; i++ )
346  {
347  currentCache->ringCounts[i] = mpoly->geoms[i]->nrings;
348  nrings += mpoly->geoms[i]->nrings;
349  }
350  currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings);
351  /*
352  ** Load the array in geometry order, each outer ring followed by the inner rings
353  ** associated with that outer ring
354  */
355  i = 0;
356  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  i++;
362  }
363  }
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  currentCache->ringCounts = lwalloc(sizeof(int));
373  currentCache->ringCounts[0] = poly->nrings;
374  /*
375  ** Just load the rings on in order
376  */
377  currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
378  for ( i = 0; i < poly->nrings; i++ )
379  {
380  currentCache->ringIndices[i] = RTreeCreate(poly->rings[i]);
381  }
382  rtree_cache->index = currentCache;
383  }
384  else
385  {
386  /* Uh oh, shouldn't be here. */
387  lwpgerror("RTreeBuilder got asked to build index on non-polygon");
388  return LW_FAILURE;
389  }
390  return LW_SUCCESS;
391 }
392 
397 static int
398 RTreeFreer(GeomCache* cache)
399 {
400  RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
401 
402  if ( ! cache )
403  return LW_FAILURE;
404 
405  if ( rtree_cache->index )
406  {
407  RTreeCacheClear(rtree_cache->index);
408  lwfree(rtree_cache->index);
409  rtree_cache->index = 0;
410  rtree_cache->argnum = 0;
411  }
412  return LW_SUCCESS;
413 }
414 
415 static GeomCache*
417 {
418  RTreeGeomCache* cache = palloc(sizeof(RTreeGeomCache));
419  memset(cache, 0, sizeof(RTreeGeomCache));
420  return (GeomCache*)cache;
421 }
422 
423 static GeomCacheMethods RTreeCacheMethods =
424 {
425  RTREE_CACHE_ENTRY,
426  RTreeBuilder,
427  RTreeFreer,
429 };
430 
432 GetRtreeCache(FunctionCallInfoData* fcinfo, GSERIALIZED* g1)
433 {
434  RTreeGeomCache* cache = (RTreeGeomCache*)GetGeomCache(fcinfo, &RTreeCacheMethods, g1, NULL);
435  RTREE_POLY_CACHE* index = NULL;
436 
437  if ( cache )
438  index = cache->index;
439 
440  return index;
441 }
442 
443 
451 {
452  LWMLINE *tmp, *result;
453  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  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: not contained.", root);
462 
463  return NULL;
464  }
465 
466  /* If there is a segment defined for this node, include it. */
467  if (root->segment)
468  {
469  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: adding segment %p %d.", root, root->segment, root->segment->type);
470 
471  lwgeoms = lwalloc(sizeof(LWGEOM *));
472  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  /* If there is a left child node, recursively include its results. */
480  if (root->leftNode)
481  {
482  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing left.", root);
483 
484  tmp = RTreeFindLineSegments(root->leftNode, value);
485  if (tmp)
486  {
487  POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));
488 
489  if (result)
490  result = RTreeMergeMultiLines(result, tmp);
491  else
492  result = tmp;
493  }
494  }
495 
496  /* Same for any right child. */
497  if (root->rightNode)
498  {
499  POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing right.", root);
500 
501  tmp = RTreeFindLineSegments(root->rightNode, value);
502  if (tmp)
503  {
504  POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, FLAGS_GET_Z(tmp->flags));
505 
506  if (result)
507  result = RTreeMergeMultiLines(result, tmp);
508  else
509  result = tmp;
510  }
511  }
512 
513  return result;
514 }
515 
516 
517 
static int RTreeFreer(GeomCache *cache)
Callback function sent into the GetGeomCache generic caching system.
Definition: lwgeom_rtree.c:398
uint8_t type
Definition: liblwgeom.h:417
uint8_t type
Definition: liblwgeom.h:476
static RTREE_NODE * RTreeCreateLeafNode(POINTARRAY *pa, uint32_t startPoint)
Creates a leaf node given the pointer to the start point of the segment.
Definition: lwgeom_rtree.c:169
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:43
char * r
Definition: cu_in_wkt.c:24
static LWMLINE * RTreeMergeMultiLines(LWMLINE *line1, LWMLINE *line2)
Merges two multilinestrings into a single multilinestring.
Definition: lwgeom_rtree.c:282
void lwfree(void *mem)
Definition: lwutil.c:244
uint8_t type
Definition: liblwgeom.h:502
uint32_t ngeoms
Definition: liblwgeom.h:493
static GeomCache * RTreeAllocator(void)
Definition: lwgeom_rtree.c:416
The following struct and methods are used for a 1D RTree implementation, described at: http://lin-ear...
Definition: lwgeom_rtree.h:45
#define POLYGONTYPE
Definition: liblwgeom.h:86
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
The tree structure used for fast P-i-P tests by point_in_multipolygon_rtree()
Definition: lwgeom_rtree.h:57
#define LW_SUCCESS
Definition: liblwgeom.h:79
void lwline_free(LWLINE *line)
Definition: lwline.c:76
struct rtree_node * leftNode
Definition: lwgeom_rtree.h:48
LWMLINE * RTreeFindLineSegments(RTREE_NODE *root, double value)
Retrieves a collection of line segments given the root and crossing value.
Definition: lwgeom_rtree.c:450
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:56
#define FP_MIN(A, B)
Representation for the y-axis interval spanned by an edge.
Definition: lwgeom_rtree.h:33
uint32_t ngeoms
Definition: liblwgeom.h:506
#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:101
uint32_t nrings
Definition: liblwgeom.h:454
#define LW_FAILURE
Definition: liblwgeom.h:78
static RTREE_NODE * RTreeCreateInteriorNode(RTREE_NODE *left, RTREE_NODE *right)
Creates an interior node given the children.
Definition: lwgeom_rtree.c:148
RTREE_INTERVAL * interval
Definition: lwgeom_rtree.h:47
unsigned int uint32_t
Definition: uthash.h:78
uint8_t flags
Definition: liblwgeom.h:477
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:110
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_FALSE, then a duplicate point will not be added.
Definition: ptarray.c:156
static RTREE_NODE * RTreeCreate(POINTARRAY *pointArray)
Creates an rtree given a pointer to the point array.
Definition: lwgeom_rtree.c:218
LWPOLY ** geoms
Definition: liblwgeom.h:495
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
static RTREE_POLY_CACHE * RTreeCacheCreate()
Allocate a fresh clean RTREE_POLY_CACHE.
Definition: lwgeom_rtree.c:43
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
LWGEOM ** lwgeoms
Definition: geobuf.h:53
POINTARRAY ** rings
Definition: liblwgeom.h:456
static int RTreeBuilder(const LWGEOM *lwgeom, GeomCache *cache)
Callback function sent into the GetGeomCache generic caching system.
Definition: lwgeom_rtree.c:316
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:139
static GeomCacheMethods RTreeCacheMethods
Definition: lwgeom_rtree.c:423
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:475
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
uint32_t ngeoms
Definition: liblwgeom.h:480
LWLINE * segment
Definition: lwgeom_rtree.h:50
LWLINE ** geoms
Definition: liblwgeom.h:482
static void RTreeCacheClear(RTREE_POLY_CACHE *cache)
Free the cache object and all the sub-objects properly.
Definition: lwgeom_rtree.c:76
struct rtree_node * rightNode
Definition: lwgeom_rtree.h:49
static RTREE_INTERVAL * RTreeCreateInterval(double value1, double value2)
Creates an interval given the min and max values, in arbitrary order.
Definition: lwgeom_rtree.c:129
uint8_t type
Definition: liblwgeom.h:395
void * lwalloc(size_t size)
Definition: lwutil.c:229
RTREE_NODE ** ringIndices
Definition: lwgeom_rtree.h:59
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:432
double y
Definition: liblwgeom.h:351
#define MULTILINETYPE
Definition: liblwgeom.h:88
uint8_t flags
Definition: liblwgeom.h:418
RTREE_POLY_CACHE * index
Definition: lwgeom_rtree.h:74
This library is the generic geometry handling section of PostGIS.
#define FP_MAX(A, B)
uint32_t npoints
Definition: liblwgeom.h:370