PostGIS  2.2.7dev-r@@SVN_REVISION@@
lwgeom_geos_cluster.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright 2015 Daniel Baston <dbaston@gmail.com>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include <string.h>
14 #include "liblwgeom.h"
15 #include "liblwgeom_internal.h"
16 #include "lwgeom_log.h"
17 #include "lwgeom_geos.h"
18 #include "lwunionfind.h"
19 
20 static const int STRTREE_NODE_CAPACITY = 10;
21 
22 /* Utility struct used to pass information to the GEOSSTRtree_query callback */
24 {
26  char error;
27  uint32_t* p;
28  const GEOSPreparedGeometry* prep;
29  GEOSGeometry** geoms;
30 };
31 
32 /* Utility struct used to pass information to the GEOSSTRtree_query callback */
34 {
36  char error;
37  uint32_t* p;
39  double tolerance;
40 };
41 
42 /* Utility struct to keep GEOSSTRtree and associated structures to be freed after use */
43 struct STRTree
44 {
45  GEOSSTRtree* tree;
46  GEOSGeometry** envelopes;
47  uint32_t* geom_ids;
48  uint32_t num_geoms;
49 };
50 
51 static struct STRTree make_strtree(void** geoms, uint32_t num_geoms, char is_lwgeom);
52 static void destroy_strtree(struct STRTree tree);
53 static void union_if_intersecting(void* item, void* userdata);
54 static void union_if_dwithin(void* item, void* userdata);
55 static int union_intersecting_pairs(GEOSGeometry** geoms, uint32_t num_geoms, UNIONFIND* uf);
56 static int combine_geometries(UNIONFIND* uf, void** geoms, uint32_t num_geoms, void*** clustersGeoms, uint32_t* num_clusters, char is_lwgeom);
57 
59 static struct STRTree
60 make_strtree(void** geoms, uint32_t num_geoms, char is_lwgeom)
61 {
62  struct STRTree tree;
63  tree.tree = GEOSSTRtree_create(STRTREE_NODE_CAPACITY);
64  if (tree.tree == NULL)
65  {
66  return tree;
67  }
68  tree.envelopes = lwalloc(num_geoms * sizeof(GEOSGeometry*));
69  tree.geom_ids = lwalloc(num_geoms * sizeof(uint32_t));
70  tree.num_geoms = num_geoms;
71 
72  size_t i;
73  for (i = 0; i < num_geoms; i++)
74  {
75  tree.geom_ids[i] = i;
76  if (!is_lwgeom)
77  {
78  tree.envelopes[i] = GEOSEnvelope(geoms[i]);
79  }
80  else
81  {
82  const GBOX* box = lwgeom_get_bbox(geoms[i]);
83  if (box)
84  {
85  tree.envelopes[i] = GBOX2GEOS(box);
86  }
87  else
88  {
89  /* Empty geometry */
90  tree.envelopes[i] = GEOSGeom_createEmptyPolygon();
91  }
92  }
93  GEOSSTRtree_insert(tree.tree, tree.envelopes[i], &(tree.geom_ids[i]));
94  }
95  return tree;
96 }
97 
99 static void
101 {
102  size_t i;
103  GEOSSTRtree_destroy(tree.tree);
104 
105  for (i = 0; i < tree.num_geoms; i++)
106  {
107  GEOSGeom_destroy(tree.envelopes[i]);
108  }
109  lwfree(tree.geom_ids);
110  lwfree(tree.envelopes);
111 }
112 
113 /* Callback function for GEOSSTRtree_query */
114 static void
115 union_if_intersecting(void* item, void* userdata)
116 {
117  struct UnionIfIntersectingContext *cxt = userdata;
118  if (cxt->error)
119  {
120  return;
121  }
122  uint32_t q = *((uint32_t*) item);
123  uint32_t p = *(cxt->p);
124 
125  if (p != q && UF_find(cxt->uf, p) != UF_find(cxt->uf, q))
126  {
127  int geos_type = GEOSGeomTypeId(cxt->geoms[p]);
128  int geos_result;
129 
130  /* Don't build prepared a geometry around a Point or MultiPoint -
131  * there are some problems in the implementation, and it's not clear
132  * there would be a performance benefit in any case. (See #3433)
133  */
134  if (geos_type != GEOS_POINT && geos_type != GEOS_MULTIPOINT)
135  {
136  /* Lazy initialize prepared geometry */
137  if (cxt->prep == NULL)
138  {
139  cxt->prep = GEOSPrepare(cxt->geoms[p]);
140  }
141  geos_result = GEOSPreparedIntersects(cxt->prep, cxt->geoms[q]);
142  }
143  else
144  {
145  geos_result = GEOSIntersects(cxt->geoms[p], cxt->geoms[q]);
146  }
147  if (geos_result > 1)
148  {
149  cxt->error = geos_result;
150  return;
151  }
152  if (geos_result)
153  {
154  UF_union(cxt->uf, p, q);
155  }
156  }
157 }
158 
159 /* Callback function for GEOSSTRtree_query */
160 static void
161 union_if_dwithin(void* item, void* userdata)
162 {
163  struct UnionIfDWithinContext *cxt = userdata;
164  if (cxt->error)
165  {
166  return;
167  }
168  uint32_t q = *((uint32_t*) item);
169  uint32_t p = *(cxt->p);
170 
171  if (p != q && UF_find(cxt->uf, p) != UF_find(cxt->uf, q))
172  {
173  double mindist = lwgeom_mindistance2d_tolerance(cxt->geoms[p], cxt->geoms[q], cxt->tolerance);
174  if (mindist == FLT_MAX)
175  {
176  cxt->error = 1;
177  return;
178  }
179 
180  if (mindist <= cxt->tolerance)
181  {
182  UF_union(cxt->uf, p, q);
183  }
184  }
185 }
186 
187 /* Identify intersecting geometries and mark them as being in the same set */
188 static int
189 union_intersecting_pairs(GEOSGeometry** geoms, uint32_t num_geoms, UNIONFIND* uf)
190 {
191  uint32_t i;
192 
193  if (num_geoms <= 1)
194  {
195  return LW_SUCCESS;
196  }
197 
198  struct STRTree tree = make_strtree((void**) geoms, num_geoms, 0);
199  if (tree.tree == NULL)
200  {
201  destroy_strtree(tree);
202  return LW_FAILURE;
203  }
204  for (i = 0; i < num_geoms; i++)
205  {
206  if (GEOSisEmpty(geoms[i]))
207  {
208  continue;
209  }
210 
211  struct UnionIfIntersectingContext cxt =
212  {
213  .uf = uf,
214  .error = 0,
215  .p = &i,
216  .prep = NULL,
217  .geoms = geoms
218  };
219  GEOSGeometry* query_envelope = GEOSEnvelope(geoms[i]);
220  GEOSSTRtree_query(tree.tree, query_envelope, &union_if_intersecting, &cxt);
221 
222  GEOSGeom_destroy(query_envelope);
223  GEOSPreparedGeom_destroy(cxt.prep);
224  if (cxt.error)
225  {
226  return LW_FAILURE;
227  }
228  }
229 
230  destroy_strtree(tree);
231  return LW_SUCCESS;
232 }
233 
234 /* Identify geometries within a distance tolerance and mark them as being in the same set */
235 static int
236 union_pairs_within_distance(LWGEOM** geoms, uint32_t num_geoms, UNIONFIND* uf, double tolerance)
237 {
238  uint32_t i;
239 
240  if (num_geoms <= 1)
241  {
242  return LW_SUCCESS;
243  }
244 
245  struct STRTree tree = make_strtree((void**) geoms, num_geoms, 1);
246  if (tree.tree == NULL)
247  {
248  destroy_strtree(tree);
249  return LW_FAILURE;
250  }
251 
252  for (i = 0; i < num_geoms; i++)
253  {
254  struct UnionIfDWithinContext cxt =
255  {
256  .uf = uf,
257  .error = 0,
258  .p = &i,
259  .geoms = geoms,
260  .tolerance = tolerance
261  };
262 
263  const GBOX* geom_extent = lwgeom_get_bbox(geoms[i]);
264  if (!geom_extent)
265  {
266  /* Empty geometry */
267  continue;
268  }
269  GBOX* query_extent = gbox_clone(geom_extent);
270  gbox_expand(query_extent, tolerance);
271  GEOSGeometry* query_envelope = GBOX2GEOS(query_extent);
272 
273  if (!query_envelope)
274  {
275  destroy_strtree(tree);
276  return LW_FAILURE;
277  }
278 
279  GEOSSTRtree_query(tree.tree, query_envelope, &union_if_dwithin, &cxt);
280 
281  lwfree(query_extent);
282  GEOSGeom_destroy(query_envelope);
283  if (cxt.error)
284  {
285  return LW_FAILURE;
286  }
287  }
288 
289  destroy_strtree(tree);
290  return LW_SUCCESS;
291 }
292 
296 int
297 cluster_intersecting(GEOSGeometry** geoms, uint32_t num_geoms, GEOSGeometry*** clusterGeoms, uint32_t* num_clusters)
298 {
299  int cluster_success;
300  UNIONFIND* uf = UF_create(num_geoms);
301 
302  if (union_intersecting_pairs(geoms, num_geoms, uf) == LW_FAILURE)
303  {
304  UF_destroy(uf);
305  return LW_FAILURE;
306  }
307 
308  cluster_success = combine_geometries(uf, (void**) geoms, num_geoms, (void***) clusterGeoms, num_clusters, 0);
309  UF_destroy(uf);
310  return cluster_success;
311 }
312 
316 int
317 cluster_within_distance(LWGEOM** geoms, uint32_t num_geoms, double tolerance, LWGEOM*** clusterGeoms, uint32_t* num_clusters)
318 {
319  int cluster_success;
320  UNIONFIND* uf = UF_create(num_geoms);
321 
322  if (union_pairs_within_distance(geoms, num_geoms, uf, tolerance) == LW_FAILURE)
323  {
324  UF_destroy(uf);
325  return LW_FAILURE;
326  }
327 
328  cluster_success = combine_geometries(uf, (void**) geoms, num_geoms, (void***) clusterGeoms, num_clusters, 1);
329  UF_destroy(uf);
330  return cluster_success;
331 }
332 
336 static int
337 combine_geometries(UNIONFIND* uf, void** geoms, uint32_t num_geoms, void*** clusterGeoms, uint32_t* num_clusters, char is_lwgeom)
338 {
339  size_t i, j, k;
340 
341  /* Combine components of each cluster into their own GeometryCollection */
342  *num_clusters = uf->num_clusters;
343  *clusterGeoms = lwalloc(*num_clusters * sizeof(void*));
344 
345  void** geoms_in_cluster = lwalloc(num_geoms * sizeof(void*));
346  uint32_t* ordered_components = UF_ordered_by_cluster(uf);
347  for (i = 0, j = 0, k = 0; i < num_geoms; i++)
348  {
349  geoms_in_cluster[j++] = geoms[ordered_components[i]];
350 
351  /* Is this the last geometry in the component? */
352  if ((i == num_geoms - 1) ||
353  (UF_find(uf, ordered_components[i]) != UF_find(uf, ordered_components[i+1])))
354  {
355  if (k >= uf->num_clusters) {
356  /* Should not get here - it means that we have more clusters than uf->num_clusters thinks we should. */
357  return LW_FAILURE;
358  }
359 
360  if (is_lwgeom)
361  {
362  LWGEOM** components = lwalloc(j * sizeof(LWGEOM*));
363  memcpy(components, geoms_in_cluster, j * sizeof(LWGEOM*));
364  (*clusterGeoms)[k++] = lwcollection_construct(COLLECTIONTYPE, components[0]->srid, NULL, j, (LWGEOM**) components);
365  }
366  else
367  {
368  (*clusterGeoms)[k++] = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, (GEOSGeometry**) geoms_in_cluster, j);
369  }
370  j = 0;
371  }
372  }
373 
374  lwfree(geoms_in_cluster);
375  lwfree(ordered_components);
376 
377  return LW_SUCCESS;
378 }
void UF_destroy(UNIONFIND *uf)
Definition: lwunionfind.c:40
GEOSGeometry ** envelopes
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
void lwfree(void *mem)
Definition: lwutil.c:214
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:199
void gbox_expand(GBOX *g, double d)
Move the box minimums down and the maximums up by the distance provided.
Definition: g_box.c:93
static void destroy_strtree(struct STRTree tree)
Clean up STRTree after use.
const GEOSPreparedGeometry * prep
#define LW_SUCCESS
Definition: liblwgeom.h:65
GEOSGeometry * GBOX2GEOS(const GBOX *box)
static void union_if_intersecting(void *item, void *userdata)
GEOSSTRtree * tree
int cluster_intersecting(GEOSGeometry **geoms, uint32_t num_geoms, GEOSGeometry ***clusterGeoms, uint32_t *num_clusters)
Takes an array of GEOSGeometry* and constructs an array of GEOSGeometry*, where each element in the c...
#define LW_FAILURE
Definition: liblwgeom.h:64
static const int STRTREE_NODE_CAPACITY
static void union_if_dwithin(void *item, void *userdata)
uint32_t UF_find(UNIONFIND *uf, uint32_t i)
Definition: lwunionfind.c:48
static struct STRTree make_strtree(void **geoms, uint32_t num_geoms, char is_lwgeom)
Make a GEOSSTRtree of either GEOSGeometry* or LWGEOM* pointers.
int cluster_within_distance(LWGEOM **geoms, uint32_t num_geoms, double tolerance, LWGEOM ***clusterGeoms, uint32_t *num_clusters)
Takes an array of LWGEOM* and constructs an array of LWGEOM*, where each element in the constructed a...
uint32_t * UF_ordered_by_cluster(UNIONFIND *uf)
Definition: lwunionfind.c:93
const GBOX * lwgeom_get_bbox(const LWGEOM *lwgeom)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:640
uint32_t * geom_ids
uint32_t num_clusters
Definition: lwunionfind.h:22
UNIONFIND * UF_create(uint32_t N)
Definition: lwunionfind.c:21
static int union_pairs_within_distance(LWGEOM **geoms, uint32_t num_geoms, UNIONFIND *uf, double tolerance)
GBOX * gbox_clone(const GBOX *gbox)
Definition: g_box.c:41
void UF_union(UNIONFIND *uf, uint32_t i, uint32_t j)
Definition: lwunionfind.c:65
void * lwalloc(size_t size)
Definition: lwutil.c:199
uint32_t num_geoms
static int combine_geometries(UNIONFIND *uf, void **geoms, uint32_t num_geoms, void ***clustersGeoms, uint32_t *num_clusters, char is_lwgeom)
Uses a UNIONFIND to identify the set with which each input geometry is associated, and groups the geometries into GeometryCollections.
static int union_intersecting_pairs(GEOSGeometry **geoms, uint32_t num_geoms, UNIONFIND *uf)
#define COLLECTIONTYPE
Definition: liblwgeom.h:76
This library is the generic geometry handling section of PostGIS.