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 2009-2014 Sandro Santilli <strk@kbt.io>
22 * Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
23 * Copyright 2001-2003 Refractions Research Inc.
24 *
25 **********************************************************************/
26
27
28 #include "../postgis_config.h"
29
30 /* PostgreSQL */
31 #include "postgres.h"
32 #include "funcapi.h"
33 #include "utils/array.h"
34 #include "utils/builtins.h"
35 #include "utils/lsyscache.h"
36 #include "utils/numeric.h"
37 #include "access/htup_details.h"
38
39 /* PostGIS */
40 #include "lwgeom_functions_analytic.h" /* for point_in_polygon */
41 #include "lwgeom_geos.h"
42 #include "liblwgeom.h"
43 #include "lwgeom_rtree.h"
44 #include "lwgeom_geos_prepared.h"
45 #include "lwgeom_accum.h"
46
47
48 /* Return NULL on GEOS error
49 *
50 * Prints error message only if it was not for interruption, in which
51 * case we let PostgreSQL deal with the error.
52 */
53 #define HANDLE_GEOS_ERROR(label) \
54 { \
55 if (strstr(lwgeom_geos_errmsg, "InterruptedException")) \
56 ereport(ERROR, \
57 (errcode(ERRCODE_QUERY_CANCELED), errmsg("canceling statement due to user request"))); \
58 else \
59 lwpgerror("%s: %s", (label), lwgeom_geos_errmsg); \
60 PG_RETURN_NULL(); \
61 }
62
63 /*
64 ** Prototypes for SQL-bound functions
65 */
66 Datum relate_full(PG_FUNCTION_ARGS);
67 Datum relate_pattern(PG_FUNCTION_ARGS);
68 Datum disjoint(PG_FUNCTION_ARGS);
69 Datum touches(PG_FUNCTION_ARGS);
70 Datum ST_Intersects(PG_FUNCTION_ARGS);
71 Datum crosses(PG_FUNCTION_ARGS);
72 Datum contains(PG_FUNCTION_ARGS);
73 Datum within(PG_FUNCTION_ARGS);
74 Datum containsproperly(PG_FUNCTION_ARGS);
75 Datum covers(PG_FUNCTION_ARGS);
76 Datum overlaps(PG_FUNCTION_ARGS);
77 Datum isvalid(PG_FUNCTION_ARGS);
78 Datum isvalidreason(PG_FUNCTION_ARGS);
79 Datum isvaliddetail(PG_FUNCTION_ARGS);
80 Datum buffer(PG_FUNCTION_ARGS);
81 Datum ST_Intersection(PG_FUNCTION_ARGS);
82 Datum convexhull(PG_FUNCTION_ARGS);
83 Datum topologypreservesimplify(PG_FUNCTION_ARGS);
84 Datum ST_Difference(PG_FUNCTION_ARGS);
85 Datum boundary(PG_FUNCTION_ARGS);
86 Datum ST_SymDifference(PG_FUNCTION_ARGS);
87 Datum ST_Union(PG_FUNCTION_ARGS);
88 Datum issimple(PG_FUNCTION_ARGS);
89 Datum isring(PG_FUNCTION_ARGS);
90 Datum pointonsurface(PG_FUNCTION_ARGS);
91 Datum GEOSnoop(PG_FUNCTION_ARGS);
92 Datum postgis_geos_version(PG_FUNCTION_ARGS);
93 Datum centroid(PG_FUNCTION_ARGS);
94 Datum polygonize_garray(PG_FUNCTION_ARGS);
95 Datum clusterintersecting_garray(PG_FUNCTION_ARGS);
96 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS);
97 Datum linemerge(PG_FUNCTION_ARGS);
98 Datum coveredby(PG_FUNCTION_ARGS);
99 Datum hausdorffdistance(PG_FUNCTION_ARGS);
100 Datum hausdorffdistancedensify(PG_FUNCTION_ARGS);
101 Datum ST_FrechetDistance(PG_FUNCTION_ARGS);
102 Datum ST_UnaryUnion(PG_FUNCTION_ARGS);
103 Datum ST_Equals(PG_FUNCTION_ARGS);
104 Datum ST_BuildArea(PG_FUNCTION_ARGS);
105 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS);
106 Datum ST_MaximumInscribedCircle(PG_FUNCTION_ARGS);
107
108 Datum pgis_union_geometry_array(PG_FUNCTION_ARGS);
109 Datum pgis_geometry_union_finalfn(PG_FUNCTION_ARGS);
110
111 /*
112 ** Prototypes end
113 */
114
115
116 PG_FUNCTION_INFO_V1(postgis_geos_version);
postgis_geos_version(PG_FUNCTION_ARGS)117 Datum postgis_geos_version(PG_FUNCTION_ARGS)
118 {
119 const char *ver = lwgeom_geos_version();
120 text *result = cstring_to_text(ver);
121 PG_RETURN_POINTER(result);
122 }
123
124 static char
is_poly(const GSERIALIZED * g)125 is_poly(const GSERIALIZED* g)
126 {
127 int type = gserialized_get_type(g);
128 return type == POLYGONTYPE || type == MULTIPOLYGONTYPE;
129 }
130
131 static char
is_point(const GSERIALIZED * g)132 is_point(const GSERIALIZED* g)
133 {
134 int type = gserialized_get_type(g);
135 return type == POINTTYPE || type == MULTIPOINTTYPE;
136 }
137
138 /* Utility function that checks a LWPOINT and a GSERIALIZED poly against a cache.
139 * Serialized poly may be a multipart.
140 */
141 static int
pip_short_circuit(RTREE_POLY_CACHE * poly_cache,LWPOINT * point,const GSERIALIZED * gpoly)142 pip_short_circuit(RTREE_POLY_CACHE *poly_cache, LWPOINT *point, const GSERIALIZED *gpoly)
143 {
144 int result;
145
146 if ( poly_cache && poly_cache->ringIndices )
147 {
148 result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
149 }
150 else
151 {
152 LWGEOM* poly = lwgeom_from_gserialized(gpoly);
153 if ( lwgeom_get_type(poly) == POLYGONTYPE )
154 {
155 result = point_in_polygon(lwgeom_as_lwpoly(poly), point);
156 }
157 else
158 {
159 result = point_in_multipolygon(lwgeom_as_lwmpoly(poly), point);
160 }
161 lwgeom_free(poly);
162 }
163
164 return result;
165 }
166
167 /**
168 * @brief Compute the Hausdorff distance thanks to the corresponding GEOS function
169 * @example ST_HausdorffDistance {@link #hausdorffdistance} - SELECT ST_HausdorffDistance(
170 * 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
171 * 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry);
172 */
173
174 PG_FUNCTION_INFO_V1(hausdorffdistance);
hausdorffdistance(PG_FUNCTION_ARGS)175 Datum hausdorffdistance(PG_FUNCTION_ARGS)
176 {
177 GSERIALIZED *geom1;
178 GSERIALIZED *geom2;
179 GEOSGeometry *g1;
180 GEOSGeometry *g2;
181 double result;
182 int retcode;
183
184 geom1 = PG_GETARG_GSERIALIZED_P(0);
185 geom2 = PG_GETARG_GSERIALIZED_P(1);
186
187 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
188 PG_RETURN_NULL();
189
190 initGEOS(lwpgnotice, lwgeom_geos_error);
191
192 g1 = POSTGIS2GEOS(geom1);
193 if (!g1)
194 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
195
196 g2 = POSTGIS2GEOS(geom2);
197 if (!g2)
198 {
199 GEOSGeom_destroy(g1);
200 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
201 }
202
203 retcode = GEOSHausdorffDistance(g1, g2, &result);
204 GEOSGeom_destroy(g1);
205 GEOSGeom_destroy(g2);
206
207 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSHausdorffDistance");
208
209 PG_FREE_IF_COPY(geom1, 0);
210 PG_FREE_IF_COPY(geom2, 1);
211
212 PG_RETURN_FLOAT8(result);
213 }
214
215 /**
216 * @brief Compute the Hausdorff distance with densification thanks to the corresponding GEOS function
217 * @example hausdorffdistancedensify {@link #hausdorffdistancedensify} - SELECT ST_HausdorffDistance(
218 * 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
219 * 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry,
220 * 0.5);
221 */
222
223 PG_FUNCTION_INFO_V1(hausdorffdistancedensify);
hausdorffdistancedensify(PG_FUNCTION_ARGS)224 Datum hausdorffdistancedensify(PG_FUNCTION_ARGS)
225 {
226 GSERIALIZED *geom1;
227 GSERIALIZED *geom2;
228 GEOSGeometry *g1;
229 GEOSGeometry *g2;
230 double densifyFrac;
231 double result;
232 int retcode;
233
234 geom1 = PG_GETARG_GSERIALIZED_P(0);
235 geom2 = PG_GETARG_GSERIALIZED_P(1);
236 densifyFrac = PG_GETARG_FLOAT8(2);
237
238 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
239 PG_RETURN_NULL();
240
241 initGEOS(lwpgnotice, lwgeom_geos_error);
242
243 g1 = POSTGIS2GEOS(geom1);
244 if (!g1)
245 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
246
247 g2 = POSTGIS2GEOS(geom2);
248 if (!g2)
249 {
250 GEOSGeom_destroy(g1);
251 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
252 }
253
254 retcode = GEOSHausdorffDistanceDensify(g1, g2, densifyFrac, &result);
255 GEOSGeom_destroy(g1);
256 GEOSGeom_destroy(g2);
257
258 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSHausdorffDistanceDensify");
259
260 PG_FREE_IF_COPY(geom1, 0);
261 PG_FREE_IF_COPY(geom2, 1);
262
263 PG_RETURN_FLOAT8(result);
264 }
265
266 /**
267 * @brief Compute the Frechet distance with optional densification thanks to the corresponding GEOS function
268 * @example ST_FrechetDistance {@link #frechetdistance} - SELECT ST_FrechetDistance(
269 * 'LINESTRING (0 0, 50 200, 100 0, 150 200, 200 0)'::geometry,
270 * 'LINESTRING (0 200, 200 150, 0 100, 200 50, 0 0)'::geometry, 0.5);
271 */
272
273 PG_FUNCTION_INFO_V1(ST_FrechetDistance);
ST_FrechetDistance(PG_FUNCTION_ARGS)274 Datum ST_FrechetDistance(PG_FUNCTION_ARGS)
275 {
276 #if POSTGIS_GEOS_VERSION < 37
277
278 lwpgerror("The GEOS version this PostGIS binary "
279 "was compiled against (%d) doesn't support "
280 "'GEOSFechetDistance' function (3.7.0+ required)",
281 POSTGIS_GEOS_VERSION);
282 PG_RETURN_NULL();
283
284 #else /* POSTGIS_GEOS_VERSION >= 37 */
285 GSERIALIZED *geom1;
286 GSERIALIZED *geom2;
287 GEOSGeometry *g1;
288 GEOSGeometry *g2;
289 double densifyFrac;
290 double result;
291 int retcode;
292
293 geom1 = PG_GETARG_GSERIALIZED_P(0);
294 geom2 = PG_GETARG_GSERIALIZED_P(1);
295 densifyFrac = PG_GETARG_FLOAT8(2);
296
297 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
298 PG_RETURN_NULL();
299
300 initGEOS(lwpgnotice, lwgeom_geos_error);
301
302 g1 = POSTGIS2GEOS(geom1);
303 if (!g1)
304 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
305
306 g2 = POSTGIS2GEOS(geom2);
307 if (!g2)
308 {
309 GEOSGeom_destroy(g1);
310 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
311 }
312
313 if (densifyFrac <= 0.0)
314 {
315 retcode = GEOSFrechetDistance(g1, g2, &result);
316 }
317 else
318 {
319 retcode = GEOSFrechetDistanceDensify(g1, g2, densifyFrac, &result);
320 }
321
322 GEOSGeom_destroy(g1);
323 GEOSGeom_destroy(g2);
324
325 if (retcode == 0) HANDLE_GEOS_ERROR("GEOSFrechetDistance");
326
327 PG_FREE_IF_COPY(geom1, 0);
328 PG_FREE_IF_COPY(geom2, 1);
329
330 PG_RETURN_FLOAT8(result);
331
332 #endif /* POSTGIS_GEOS_VERSION >= 37 */
333 }
334
335
336 /**
337 * @brief Compute the Frechet distance with optional densification thanks to the corresponding GEOS function
338 * @example ST_FrechetDistance {@link #frechetdistance} - SELECT ST_FrechetDistance(
339 * 'LINESTRING (0 0, 50 200, 100 0, 150 200, 200 0)'::geometry,
340 * 'LINESTRING (0 200, 200 150, 0 100, 200 50, 0 0)'::geometry, 0.5);
341 */
342
343 PG_FUNCTION_INFO_V1(ST_MaximumInscribedCircle);
ST_MaximumInscribedCircle(PG_FUNCTION_ARGS)344 Datum ST_MaximumInscribedCircle(PG_FUNCTION_ARGS)
345 {
346 #if POSTGIS_GEOS_VERSION < 39
347
348 lwpgerror("The GEOS version this PostGIS binary "
349 "was compiled against (%d) doesn't support "
350 "'GEOSMaximumInscribedCircle' function (3.9.0+ required)",
351 POSTGIS_GEOS_VERSION);
352 PG_RETURN_NULL();
353
354 #else /* POSTGIS_GEOS_VERSION >= 39 */
355 GSERIALIZED* geom;
356 GSERIALIZED* center;
357 GSERIALIZED* nearest;
358 TupleDesc resultTupleDesc;
359 HeapTuple resultTuple;
360 Datum result;
361 Datum result_values[3];
362 bool result_is_null[3];
363 double radius = 0.0;
364 int32 srid = SRID_UNKNOWN;
365 bool is3d;
366
367 if (PG_ARGISNULL(0))
368 PG_RETURN_NULL();
369
370 geom = PG_GETARG_GSERIALIZED_P(0);
371 srid = gserialized_get_srid(geom);
372 is3d = gserialized_has_z(geom);
373
374 /* Empty geometry? Return POINT EMPTY with zero radius */
375 if (gserialized_is_empty(geom))
376 {
377 LWGEOM* lwcenter = (LWGEOM*) lwpoint_construct_empty(gserialized_get_srid(geom), LW_FALSE, LW_FALSE);
378 LWGEOM* lwnearest = (LWGEOM*) lwpoint_construct_empty(gserialized_get_srid(geom), LW_FALSE, LW_FALSE);
379 center = geometry_serialize(lwcenter);
380 nearest = geometry_serialize(lwnearest);
381 radius = 0.0;
382 }
383 else
384 {
385 GEOSGeometry *ginput, *gcircle, *gcenter, *gnearest;
386 double width, height, size, tolerance;
387 GBOX gbox;
388 int gtype;
389
390 if (!gserialized_get_gbox_p(geom, &gbox))
391 PG_RETURN_NULL();
392
393 width = gbox.xmax - gbox.xmin;
394 height = gbox.ymax - gbox.ymin;
395 size = width > height ? width : height;
396 tolerance = size / 1000.0;
397
398 initGEOS(lwpgnotice, lwgeom_geos_error);
399
400 ginput = POSTGIS2GEOS(geom);
401 if (!ginput)
402 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
403
404 gtype = gserialized_get_type(geom);
405 if (gtype == POLYGONTYPE || gtype == MULTIPOLYGONTYPE)
406 {
407 gcircle = GEOSMaximumInscribedCircle(ginput, tolerance);
408 if (!gcircle)
409 {
410 lwpgerror("Error calculating GEOSMaximumInscribedCircle.");
411 GEOSGeom_destroy(ginput);
412 PG_RETURN_NULL();
413 }
414 }
415 else
416 {
417 gcircle = GEOSLargestEmptyCircle(ginput, NULL, tolerance);
418 if (!gcircle)
419 {
420 lwpgerror("Error calculating GEOSLargestEmptyCircle.");
421 GEOSGeom_destroy(ginput);
422 PG_RETURN_NULL();
423 }
424 }
425
426 gcenter = GEOSGeomGetStartPoint(gcircle);
427 gnearest = GEOSGeomGetEndPoint(gcircle);
428 GEOSDistance(gcenter, gnearest, &radius);
429 GEOSSetSRID(gcenter, srid);
430 GEOSSetSRID(gnearest, srid);
431
432 center = GEOS2POSTGIS(gcenter, is3d);
433 nearest = GEOS2POSTGIS(gnearest, is3d);
434 GEOSGeom_destroy(gcenter);
435 GEOSGeom_destroy(gnearest);
436 GEOSGeom_destroy(gcircle);
437 GEOSGeom_destroy(ginput);
438 }
439
440 get_call_result_type(fcinfo, NULL, &resultTupleDesc);
441 BlessTupleDesc(resultTupleDesc);
442
443 result_values[0] = PointerGetDatum(center);
444 result_is_null[0] = false;
445 result_values[1] = PointerGetDatum(nearest);
446 result_is_null[1] = false;
447 result_values[2] = Float8GetDatum(radius);
448 result_is_null[2] = false;
449 resultTuple = heap_form_tuple(resultTupleDesc, result_values, result_is_null);
450
451 result = HeapTupleGetDatum(resultTuple);
452
453 PG_RETURN_DATUM(result);
454
455 #endif /* POSTGIS_GEOS_VERSION >= 39 */
456 }
457
458
459 /**
460 * @brief This is the final function for GeomUnion
461 * aggregate. Will have as input an array of Geometries.
462 * Will iteratively call GEOSUnion on the GEOS-converted
463 * versions of them and return PGIS-converted version back.
464 * Changing combination order *might* speed up performance.
465 */
466 PG_FUNCTION_INFO_V1(pgis_union_geometry_array);
pgis_union_geometry_array(PG_FUNCTION_ARGS)467 Datum pgis_union_geometry_array(PG_FUNCTION_ARGS)
468 {
469 ArrayType *array;
470
471 ArrayIterator iterator;
472 Datum value;
473 bool isnull;
474
475 int is3d = LW_FALSE, gotsrid = LW_FALSE;
476 int nelems = 0, geoms_size = 0, curgeom = 0, count = 0;
477
478 GSERIALIZED *gser_out = NULL;
479
480 GEOSGeometry *g = NULL;
481 GEOSGeometry *g_union = NULL;
482 GEOSGeometry **geoms = NULL;
483
484 int32_t srid = SRID_UNKNOWN;
485
486 int empty_type = 0;
487
488 /* Null array, null geometry (should be empty?) */
489 if ( PG_ARGISNULL(0) )
490 PG_RETURN_NULL();
491
492 array = PG_GETARG_ARRAYTYPE_P(0);
493 nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
494
495 /* Empty array? Null return */
496 if ( nelems == 0 ) PG_RETURN_NULL();
497
498 /* Quick scan for nulls */
499 iterator = array_create_iterator(array, 0, NULL);
500 while (array_iterate(iterator, &value, &isnull))
501 {
502 /* Skip null array items */
503 if (isnull) continue;
504 count++;
505 }
506 array_free_iterator(iterator);
507
508
509 /* All-nulls? Return null */
510 if ( count == 0 )
511 PG_RETURN_NULL();
512
513 /* One geom, good geom? Return it */
514 if ( count == 1 && nelems == 1 )
515 {
516 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
517 #endif
518 PG_RETURN_POINTER((GSERIALIZED *)(ARR_DATA_PTR(array)));
519 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
520 #endif
521 }
522
523 /* Ok, we really need GEOS now ;) */
524 initGEOS(lwpgnotice, lwgeom_geos_error);
525
526 /*
527 ** Collect the non-empty inputs and stuff them into a GEOS collection
528 */
529 geoms_size = nelems;
530 geoms = palloc(sizeof(GEOSGeometry*) * geoms_size);
531
532 /*
533 ** We need to convert the array of GSERIALIZED into a GEOS collection.
534 ** First make an array of GEOS geometries.
535 */
536 iterator = array_create_iterator(array, 0, NULL);
537 while (array_iterate(iterator, &value, &isnull))
538 {
539 GSERIALIZED *gser_in;
540
541 /* Skip null array items */
542 if (isnull) continue;
543 gser_in = (GSERIALIZED *)DatumGetPointer(value);
544
545 /* Check for SRID mismatch in array elements */
546 if ( gotsrid )
547 gserialized_error_if_srid_mismatch_reference(gser_in, srid, __func__);
548 else
549 {
550 /* Initialize SRID/dimensions info */
551 srid = gserialized_get_srid(gser_in);
552 is3d = gserialized_has_z(gser_in);
553 gotsrid = 1;
554 }
555
556 /* Don't include empties in the union */
557 if ( gserialized_is_empty(gser_in) )
558 {
559 int gser_type = gserialized_get_type(gser_in);
560 if (gser_type > empty_type)
561 {
562 empty_type = gser_type;
563 POSTGIS_DEBUGF(4, "empty_type = %d gser_type = %d", empty_type, gser_type);
564 }
565 }
566 else
567 {
568 g = POSTGIS2GEOS(gser_in);
569
570 /* Uh oh! Exception thrown at construction... */
571 if ( ! g )
572 {
573 HANDLE_GEOS_ERROR(
574 "One of the geometries in the set "
575 "could not be converted to GEOS");
576 }
577
578 /* Ensure we have enough space in our storage array */
579 if ( curgeom == geoms_size )
580 {
581 geoms_size *= 2;
582 geoms = repalloc( geoms, sizeof(GEOSGeometry*) * geoms_size );
583 }
584
585 geoms[curgeom] = g;
586 curgeom++;
587 }
588
589 }
590 array_free_iterator(iterator);
591
592 /*
593 ** Take our GEOS geometries and turn them into a GEOS collection,
594 ** then pass that into cascaded union.
595 */
596 if (curgeom > 0)
597 {
598 g = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, curgeom);
599 if (!g) HANDLE_GEOS_ERROR("Could not create GEOS COLLECTION from geometry array");
600
601 g_union = GEOSUnaryUnion(g);
602 GEOSGeom_destroy(g);
603 if (!g_union) HANDLE_GEOS_ERROR("GEOSUnaryUnion");
604
605 GEOSSetSRID(g_union, srid);
606 gser_out = GEOS2POSTGIS(g_union, is3d);
607 GEOSGeom_destroy(g_union);
608 }
609 /* No real geometries in our array, any empties? */
610 else
611 {
612 /* If it was only empties, we'll return the largest type number */
613 if ( empty_type > 0 )
614 {
615 PG_RETURN_POINTER(geometry_serialize(lwgeom_construct_empty(empty_type, srid, is3d, 0)));
616 }
617 /* Nothing but NULL, returns NULL */
618 else
619 {
620 PG_RETURN_NULL();
621 }
622 }
623
624 if ( ! gser_out )
625 {
626 /* Union returned a NULL geometry */
627 PG_RETURN_NULL();
628 }
629
630 PG_RETURN_POINTER(gser_out);
631 }
632
633
634 PG_FUNCTION_INFO_V1(pgis_geometry_union_finalfn);
pgis_geometry_union_finalfn(PG_FUNCTION_ARGS)635 Datum pgis_geometry_union_finalfn(PG_FUNCTION_ARGS)
636 {
637 CollectionBuildState *state;
638 ListCell *l;
639 LWGEOM **geoms;
640 GSERIALIZED *gser_out;
641 size_t ngeoms = 0;
642 int empty_type = 0;
643 bool first = true;
644 int32_t srid = SRID_UNKNOWN;
645 int has_z = LW_FALSE;
646
647 if (PG_ARGISNULL(0))
648 PG_RETURN_NULL(); /* returns null iff no input values */
649
650 state = (CollectionBuildState *)PG_GETARG_POINTER(0);
651 geoms = palloc(list_length(state->geoms) * sizeof(LWGEOM*));
652
653 /* Read contents of list into an array of only non-null values */
654 foreach (l, state->geoms)
655 {
656 LWGEOM *geom = (LWGEOM*)(lfirst(l));
657 if (geom)
658 {
659 if (!lwgeom_is_empty(geom))
660 {
661 geoms[ngeoms++] = geom;
662 if (first)
663 {
664 srid = lwgeom_get_srid(geom);
665 has_z = lwgeom_has_z(geom);
666 first = false;
667 }
668 }
669 else
670 {
671 int type = lwgeom_get_type(geom);
672 empty_type = type > empty_type ? type : empty_type;
673 srid = (srid != SRID_UNKNOWN ? srid : lwgeom_get_srid(geom));
674 }
675 }
676 }
677
678 /*
679 ** Take our array of LWGEOM* and turn it into a GEOS collection,
680 ** then pass that into cascaded union.
681 */
682 if (ngeoms > 0)
683 {
684 LWCOLLECTION* col = lwcollection_construct(COLLECTIONTYPE, srid, NULL, ngeoms, geoms);
685 LWGEOM *out = lwgeom_unaryunion_prec(lwcollection_as_lwgeom(col), state->gridSize);
686 if ( ! out )
687 {
688 lwcollection_free(col);
689 }
690 gser_out = geometry_serialize(out);
691 }
692 /* No real geometries in our array, any empties? */
693 else
694 {
695 /* If it was only empties, we'll return the largest type number */
696 if (empty_type > 0)
697 PG_RETURN_POINTER(
698 geometry_serialize(lwgeom_construct_empty(empty_type, srid, has_z, 0)));
699
700 /* Nothing but NULL, returns NULL */
701 else
702 PG_RETURN_NULL();
703 }
704
705 if (!gser_out)
706 {
707 /* Union returned a NULL geometry */
708 PG_RETURN_NULL();
709 }
710
711 PG_RETURN_POINTER(gser_out);
712 }
713
714
715
716 /**
717 * @example ST_UnaryUnion {@link #geomunion} SELECT ST_UnaryUnion(
718 * 'POLYGON((0 0, 10 0, 0 10, 10 10, 0 0))'
719 * );
720 *
721 */
722 PG_FUNCTION_INFO_V1(ST_UnaryUnion);
ST_UnaryUnion(PG_FUNCTION_ARGS)723 Datum ST_UnaryUnion(PG_FUNCTION_ARGS)
724 {
725 GSERIALIZED *geom1;
726 GSERIALIZED *result;
727 LWGEOM *lwgeom1, *lwresult ;
728 double prec = -1;
729
730 geom1 = PG_GETARG_GSERIALIZED_P(0);
731 if (PG_NARGS() > 1 && ! PG_ARGISNULL(1))
732 prec = PG_GETARG_FLOAT8(1);
733
734 lwgeom1 = lwgeom_from_gserialized(geom1) ;
735
736 lwresult = lwgeom_unaryunion_prec(lwgeom1, prec);
737 result = geometry_serialize(lwresult) ;
738
739 lwgeom_free(lwgeom1) ;
740 lwgeom_free(lwresult) ;
741
742 PG_FREE_IF_COPY(geom1, 0);
743
744 PG_RETURN_POINTER(result);
745 }
746
747 PG_FUNCTION_INFO_V1(ST_Union);
ST_Union(PG_FUNCTION_ARGS)748 Datum ST_Union(PG_FUNCTION_ARGS)
749 {
750 GSERIALIZED *geom1;
751 GSERIALIZED *geom2;
752 GSERIALIZED *result;
753 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
754 double gridSize = -1;
755
756 geom1 = PG_GETARG_GSERIALIZED_P(0);
757 geom2 = PG_GETARG_GSERIALIZED_P(1);
758 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
759 gridSize = PG_GETARG_FLOAT8(2);
760
761 lwgeom1 = lwgeom_from_gserialized(geom1);
762 lwgeom2 = lwgeom_from_gserialized(geom2);
763
764 lwresult = lwgeom_union_prec(lwgeom1, lwgeom2, gridSize);
765 result = geometry_serialize(lwresult);
766
767 lwgeom_free(lwgeom1);
768 lwgeom_free(lwgeom2);
769 lwgeom_free(lwresult);
770
771 PG_FREE_IF_COPY(geom1, 0);
772 PG_FREE_IF_COPY(geom2, 1);
773
774 PG_RETURN_POINTER(result);
775 }
776
777 /* This is retained for backward ABI compatibility
778 * with PostGIS < 3.1.0 */
779 PG_FUNCTION_INFO_V1(symdifference);
symdifference(PG_FUNCTION_ARGS)780 Datum symdifference(PG_FUNCTION_ARGS)
781 {
782 PG_RETURN_DATUM(DirectFunctionCall2(
783 ST_SymDifference,
784 PG_GETARG_DATUM(0), PG_GETARG_DATUM(1)
785 ));
786 }
787
788 /**
789 * @example symdifference {@link #symdifference} - SELECT ST_SymDifference(
790 * 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
791 * 'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
792 */
793 PG_FUNCTION_INFO_V1(ST_SymDifference);
ST_SymDifference(PG_FUNCTION_ARGS)794 Datum ST_SymDifference(PG_FUNCTION_ARGS)
795 {
796 GSERIALIZED *geom1;
797 GSERIALIZED *geom2;
798 GSERIALIZED *result;
799 LWGEOM *lwgeom1, *lwgeom2, *lwresult ;
800 double prec = -1;
801
802 geom1 = PG_GETARG_GSERIALIZED_P(0);
803 geom2 = PG_GETARG_GSERIALIZED_P(1);
804 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
805 prec = PG_GETARG_FLOAT8(2);
806
807 lwgeom1 = lwgeom_from_gserialized(geom1) ;
808 lwgeom2 = lwgeom_from_gserialized(geom2) ;
809
810 lwresult = lwgeom_symdifference_prec(lwgeom1, lwgeom2, prec);
811 result = geometry_serialize(lwresult) ;
812
813 lwgeom_free(lwgeom1) ;
814 lwgeom_free(lwgeom2) ;
815 lwgeom_free(lwresult) ;
816
817 PG_FREE_IF_COPY(geom1, 0);
818 PG_FREE_IF_COPY(geom2, 1);
819
820 PG_RETURN_POINTER(result);
821 }
822
823
824 PG_FUNCTION_INFO_V1(boundary);
boundary(PG_FUNCTION_ARGS)825 Datum boundary(PG_FUNCTION_ARGS)
826 {
827 GSERIALIZED *geom1;
828 GEOSGeometry *g1, *g3;
829 GSERIALIZED *result;
830 LWGEOM *lwgeom;
831 int32_t srid;
832
833 geom1 = PG_GETARG_GSERIALIZED_P(0);
834
835 /* Empty.Boundary() == Empty */
836 if ( gserialized_is_empty(geom1) )
837 PG_RETURN_POINTER(geom1);
838
839 srid = gserialized_get_srid(geom1);
840
841 lwgeom = lwgeom_from_gserialized(geom1);
842 if ( ! lwgeom ) {
843 lwpgerror("POSTGIS2GEOS: unable to deserialize input");
844 PG_RETURN_NULL();
845 }
846
847 /* GEOS doesn't do triangle type, so we special case that here */
848 if (lwgeom->type == TRIANGLETYPE)
849 {
850 lwgeom->type = LINETYPE;
851 result = geometry_serialize(lwgeom);
852 lwgeom_free(lwgeom);
853 PG_RETURN_POINTER(result);
854 }
855
856 initGEOS(lwpgnotice, lwgeom_geos_error);
857
858 g1 = LWGEOM2GEOS(lwgeom, 0);
859 lwgeom_free(lwgeom);
860
861 if (!g1)
862 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
863
864 g3 = GEOSBoundary(g1);
865
866 if (!g3)
867 {
868 GEOSGeom_destroy(g1);
869 HANDLE_GEOS_ERROR("GEOSBoundary");
870 }
871
872 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
873
874 GEOSSetSRID(g3, srid);
875
876 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
877
878 if (!result)
879 {
880 GEOSGeom_destroy(g1);
881 GEOSGeom_destroy(g3);
882 elog(NOTICE,
883 "GEOS2POSTGIS threw an error (result postgis geometry "
884 "formation)!");
885 PG_RETURN_NULL(); /* never get here */
886 }
887
888 GEOSGeom_destroy(g1);
889 GEOSGeom_destroy(g3);
890
891 PG_FREE_IF_COPY(geom1, 0);
892
893 PG_RETURN_POINTER(result);
894 }
895
896 PG_FUNCTION_INFO_V1(convexhull);
convexhull(PG_FUNCTION_ARGS)897 Datum convexhull(PG_FUNCTION_ARGS)
898 {
899 GSERIALIZED *geom1;
900 GEOSGeometry *g1, *g3;
901 GSERIALIZED *result;
902 LWGEOM *lwout;
903 int32_t srid;
904 GBOX bbox;
905
906 geom1 = PG_GETARG_GSERIALIZED_P(0);
907
908 /* Empty.ConvexHull() == Empty */
909 if ( gserialized_is_empty(geom1) )
910 PG_RETURN_POINTER(geom1);
911
912 srid = gserialized_get_srid(geom1);
913
914 initGEOS(lwpgnotice, lwgeom_geos_error);
915
916 g1 = POSTGIS2GEOS(geom1);
917
918 if (!g1)
919 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
920
921 g3 = GEOSConvexHull(g1);
922 GEOSGeom_destroy(g1);
923
924 if (!g3) HANDLE_GEOS_ERROR("GEOSConvexHull");
925
926 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
927
928 GEOSSetSRID(g3, srid);
929
930 lwout = GEOS2LWGEOM(g3, gserialized_has_z(geom1));
931 GEOSGeom_destroy(g3);
932
933 if (!lwout)
934 {
935 elog(ERROR,
936 "convexhull() failed to convert GEOS geometry to LWGEOM");
937 PG_RETURN_NULL(); /* never get here */
938 }
939
940 /* Copy input bbox if any */
941 if ( gserialized_get_gbox_p(geom1, &bbox) )
942 {
943 /* Force the box to have the same dimensionality as the lwgeom */
944 bbox.flags = lwout->flags;
945 lwout->bbox = gbox_copy(&bbox);
946 }
947
948 result = geometry_serialize(lwout);
949 lwgeom_free(lwout);
950
951 if (!result)
952 {
953 elog(ERROR,"GEOS convexhull() threw an error (result postgis geometry formation)!");
954 PG_RETURN_NULL(); /* never get here */
955 }
956
957 PG_FREE_IF_COPY(geom1, 0);
958 PG_RETURN_POINTER(result);
959 }
960
961 PG_FUNCTION_INFO_V1(topologypreservesimplify);
topologypreservesimplify(PG_FUNCTION_ARGS)962 Datum topologypreservesimplify(PG_FUNCTION_ARGS)
963 {
964 GSERIALIZED *geom1;
965 double tolerance;
966 GEOSGeometry *g1, *g3;
967 GSERIALIZED *result;
968 uint32_t type;
969
970 geom1 = PG_GETARG_GSERIALIZED_P(0);
971 tolerance = PG_GETARG_FLOAT8(1);
972
973 /* Empty.Simplify() == Empty */
974 type = gserialized_get_type(geom1);
975 if ( gserialized_is_empty(geom1) || type == TINTYPE || type == TRIANGLETYPE )
976 PG_RETURN_POINTER(geom1);
977
978 initGEOS(lwpgnotice, lwgeom_geos_error);
979
980 g1 = POSTGIS2GEOS(geom1);
981 if (!g1)
982 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
983
984 g3 = GEOSTopologyPreserveSimplify(g1,tolerance);
985 GEOSGeom_destroy(g1);
986
987 if (!g3) HANDLE_GEOS_ERROR("GEOSTopologyPreserveSimplify");
988
989 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
990
991 GEOSSetSRID(g3, gserialized_get_srid(geom1));
992
993 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
994 GEOSGeom_destroy(g3);
995
996 if (!result)
997 {
998 elog(ERROR,"GEOS topologypreservesimplify() threw an error (result postgis geometry formation)!");
999 PG_RETURN_NULL(); /* never get here */
1000 }
1001
1002 PG_FREE_IF_COPY(geom1, 0);
1003 PG_RETURN_POINTER(result);
1004 }
1005
1006 PG_FUNCTION_INFO_V1(buffer);
buffer(PG_FUNCTION_ARGS)1007 Datum buffer(PG_FUNCTION_ARGS)
1008 {
1009 GEOSBufferParams *bufferparams;
1010 GEOSGeometry *g1, *g3 = NULL;
1011 GSERIALIZED *result;
1012 LWGEOM *lwg;
1013 int quadsegs = 8; /* the default */
1014 int singleside = 0; /* the default */
1015 enum
1016 {
1017 ENDCAP_ROUND = 1,
1018 ENDCAP_FLAT = 2,
1019 ENDCAP_SQUARE = 3
1020 };
1021 enum
1022 {
1023 JOIN_ROUND = 1,
1024 JOIN_MITRE = 2,
1025 JOIN_BEVEL = 3
1026 };
1027 const double DEFAULT_MITRE_LIMIT = 5.0;
1028 const int DEFAULT_ENDCAP_STYLE = ENDCAP_ROUND;
1029 const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1030 double mitreLimit = DEFAULT_MITRE_LIMIT;
1031 int endCapStyle = DEFAULT_ENDCAP_STYLE;
1032 int joinStyle = DEFAULT_JOIN_STYLE;
1033
1034 GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
1035 double size = PG_GETARG_FLOAT8(1);
1036 text *params_text;
1037
1038 if (PG_NARGS() > 2)
1039 {
1040 params_text = PG_GETARG_TEXT_P(2);
1041 }
1042 else
1043 {
1044 params_text = cstring_to_text("");
1045 }
1046
1047 /* Empty.Buffer() == Empty[polygon] */
1048 if ( gserialized_is_empty(geom1) )
1049 {
1050 lwg = lwpoly_as_lwgeom(lwpoly_construct_empty(
1051 gserialized_get_srid(geom1),
1052 0, 0)); // buffer wouldn't give back z or m anyway
1053 PG_RETURN_POINTER(geometry_serialize(lwg));
1054 }
1055
1056 initGEOS(lwpgnotice, lwgeom_geos_error);
1057
1058 g1 = POSTGIS2GEOS(geom1);
1059 if (!g1)
1060 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1061
1062
1063 if (VARSIZE_ANY_EXHDR(params_text) > 0)
1064 {
1065 char *param;
1066 char *params = text_to_cstring(params_text);
1067
1068 for (param=params; ; param=NULL)
1069 {
1070 char *key, *val;
1071 param = strtok(param, " ");
1072 if (!param) break;
1073 POSTGIS_DEBUGF(3, "Param: %s", param);
1074
1075 key = param;
1076 val = strchr(key, '=');
1077 if (!val || *(val + 1) == '\0')
1078 {
1079 lwpgerror("Missing value for buffer parameter %s", key);
1080 break;
1081 }
1082 *val = '\0';
1083 ++val;
1084
1085 POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1086
1087 if ( !strcmp(key, "endcap") )
1088 {
1089 /* Supported end cap styles:
1090 * "round", "flat", "square"
1091 */
1092 if ( !strcmp(val, "round") )
1093 {
1094 endCapStyle = ENDCAP_ROUND;
1095 }
1096 else if ( !strcmp(val, "flat") ||
1097 !strcmp(val, "butt") )
1098 {
1099 endCapStyle = ENDCAP_FLAT;
1100 }
1101 else if ( !strcmp(val, "square") )
1102 {
1103 endCapStyle = ENDCAP_SQUARE;
1104 }
1105 else
1106 {
1107 lwpgerror("Invalid buffer end cap "
1108 "style: %s (accept: "
1109 "'round', 'flat', 'butt' "
1110 "or 'square'"
1111 ")", val);
1112 break;
1113 }
1114
1115 }
1116 else if ( !strcmp(key, "join") )
1117 {
1118 if ( !strcmp(val, "round") )
1119 {
1120 joinStyle = JOIN_ROUND;
1121 }
1122 else if ( !strcmp(val, "mitre") ||
1123 !strcmp(val, "miter") )
1124 {
1125 joinStyle = JOIN_MITRE;
1126 }
1127 else if ( !strcmp(val, "bevel") )
1128 {
1129 joinStyle = JOIN_BEVEL;
1130 }
1131 else
1132 {
1133 lwpgerror("Invalid buffer end cap "
1134 "style: %s (accept: "
1135 "'round', 'mitre', 'miter' "
1136 " or 'bevel'"
1137 ")", val);
1138 break;
1139 }
1140 }
1141 else if ( !strcmp(key, "mitre_limit") ||
1142 !strcmp(key, "miter_limit") )
1143 {
1144 /* mitreLimit is a float */
1145 mitreLimit = atof(val);
1146 }
1147 else if ( !strcmp(key, "quad_segs") )
1148 {
1149 /* quadrant segments is an int */
1150 quadsegs = atoi(val);
1151 }
1152 else if ( !strcmp(key, "side") )
1153 {
1154 if ( !strcmp(val, "both") )
1155 {
1156 singleside = 0;
1157 }
1158 else if ( !strcmp(val, "left") )
1159 {
1160 singleside = 1;
1161 }
1162 else if ( !strcmp(val, "right") )
1163 {
1164 singleside = 1;
1165 size *= -1;
1166 }
1167 else
1168 {
1169 lwpgerror("Invalid side parameter: %s (accept: 'right', 'left', 'both')", val);
1170 break;
1171 }
1172 }
1173 else
1174 {
1175 lwpgerror(
1176 "Invalid buffer parameter: %s (accept: 'endcap', 'join', 'mitre_limit', 'miter_limit', 'quad_segs' and 'side')",
1177 key);
1178 break;
1179 }
1180 }
1181 pfree(params); /* was pstrduped */
1182 }
1183
1184
1185 POSTGIS_DEBUGF(3, "endCap:%d joinStyle:%d mitreLimit:%g",
1186 endCapStyle, joinStyle, mitreLimit);
1187
1188 bufferparams = GEOSBufferParams_create();
1189 if (bufferparams)
1190 {
1191 if (GEOSBufferParams_setEndCapStyle(bufferparams, endCapStyle) &&
1192 GEOSBufferParams_setJoinStyle(bufferparams, joinStyle) &&
1193 GEOSBufferParams_setMitreLimit(bufferparams, mitreLimit) &&
1194 GEOSBufferParams_setQuadrantSegments(bufferparams, quadsegs) &&
1195 GEOSBufferParams_setSingleSided(bufferparams, singleside))
1196 {
1197 g3 = GEOSBufferWithParams(g1, bufferparams, size);
1198 }
1199 else
1200 {
1201 lwpgerror("Error setting buffer parameters.");
1202 }
1203 GEOSBufferParams_destroy(bufferparams);
1204 }
1205 else
1206 {
1207 lwpgerror("Error setting buffer parameters.");
1208 }
1209
1210 GEOSGeom_destroy(g1);
1211
1212 if (!g3) HANDLE_GEOS_ERROR("GEOSBuffer");
1213
1214 POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1215
1216 GEOSSetSRID(g3, gserialized_get_srid(geom1));
1217
1218 result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
1219 GEOSGeom_destroy(g3);
1220
1221 if (!result)
1222 {
1223 elog(ERROR,"GEOS buffer() threw an error (result postgis geometry formation)!");
1224 PG_RETURN_NULL(); /* never get here */
1225 }
1226
1227 PG_FREE_IF_COPY(geom1, 0);
1228 PG_RETURN_POINTER(result);
1229 }
1230
1231 /*
1232 * Generate a field of random points within the area of a
1233 * polygon or multipolygon. Throws an error for other geometry
1234 * types.
1235 */
1236 Datum ST_GeneratePoints(PG_FUNCTION_ARGS);
1237 PG_FUNCTION_INFO_V1(ST_GeneratePoints);
ST_GeneratePoints(PG_FUNCTION_ARGS)1238 Datum ST_GeneratePoints(PG_FUNCTION_ARGS)
1239 {
1240 GSERIALIZED *gser_input;
1241 GSERIALIZED *gser_result;
1242 LWGEOM *lwgeom_input;
1243 LWGEOM *lwgeom_result;
1244 int32 npoints;
1245 int32 seed = 0;
1246
1247 gser_input = PG_GETARG_GSERIALIZED_P(0);
1248 npoints = PG_GETARG_INT32(1);
1249
1250 if (npoints < 0)
1251 PG_RETURN_NULL();
1252
1253 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
1254 {
1255 seed = PG_GETARG_INT32(2);
1256 if (seed < 1)
1257 {
1258 lwpgerror("ST_GeneratePoints: seed must be greater than zero");
1259 PG_RETURN_NULL();
1260 }
1261 }
1262
1263 /* Types get checked in the code, we'll keep things small out there */
1264 lwgeom_input = lwgeom_from_gserialized(gser_input);
1265 lwgeom_result = (LWGEOM*)lwgeom_to_points(lwgeom_input, npoints, seed);
1266 lwgeom_free(lwgeom_input);
1267 PG_FREE_IF_COPY(gser_input, 0);
1268
1269 /* Return null as null */
1270 if (!lwgeom_result)
1271 PG_RETURN_NULL();
1272
1273 /* Serialize and return */
1274 gser_result = geometry_serialize(lwgeom_result);
1275 lwgeom_free(lwgeom_result);
1276 PG_RETURN_POINTER(gser_result);
1277 }
1278
1279
1280 /*
1281 * Compute at offset curve to a line
1282 */
1283 Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
1284 PG_FUNCTION_INFO_V1(ST_OffsetCurve);
ST_OffsetCurve(PG_FUNCTION_ARGS)1285 Datum ST_OffsetCurve(PG_FUNCTION_ARGS)
1286 {
1287 GSERIALIZED *gser_input;
1288 GSERIALIZED *gser_result;
1289 LWGEOM *lwgeom_input;
1290 LWGEOM *lwgeom_result;
1291 double size;
1292 int quadsegs = 8; /* the default */
1293 int nargs;
1294
1295 enum
1296 {
1297 JOIN_ROUND = 1,
1298 JOIN_MITRE = 2,
1299 JOIN_BEVEL = 3
1300 };
1301
1302 static const double DEFAULT_MITRE_LIMIT = 5.0;
1303 static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1304 double mitreLimit = DEFAULT_MITRE_LIMIT;
1305 int joinStyle = DEFAULT_JOIN_STYLE;
1306 char *param = NULL;
1307 char *paramstr = NULL;
1308
1309 /* Read SQL arguments */
1310 nargs = PG_NARGS();
1311 gser_input = PG_GETARG_GSERIALIZED_P(0);
1312 size = PG_GETARG_FLOAT8(1);
1313
1314 /* For distance == 0, just return the input. */
1315 if (size == 0) PG_RETURN_POINTER(gser_input);
1316
1317 /* Read the lwgeom, check for errors */
1318 lwgeom_input = lwgeom_from_gserialized(gser_input);
1319 if ( ! lwgeom_input )
1320 lwpgerror("ST_OffsetCurve: lwgeom_from_gserialized returned NULL");
1321
1322 /* For empty inputs, just echo them back */
1323 if ( lwgeom_is_empty(lwgeom_input) )
1324 PG_RETURN_POINTER(gser_input);
1325
1326 /* Process the optional arguments */
1327 if ( nargs > 2 )
1328 {
1329 text *wkttext = PG_GETARG_TEXT_P(2);
1330 paramstr = text_to_cstring(wkttext);
1331
1332 POSTGIS_DEBUGF(3, "paramstr: %s", paramstr);
1333
1334 for ( param=paramstr; ; param=NULL )
1335 {
1336 char *key, *val;
1337 param = strtok(param, " ");
1338 if (!param) break;
1339 POSTGIS_DEBUGF(3, "Param: %s", param);
1340
1341 key = param;
1342 val = strchr(key, '=');
1343 if (!val || *(val + 1) == '\0')
1344 {
1345 lwpgerror("ST_OffsetCurve: Missing value for buffer parameter %s", key);
1346 break;
1347 }
1348 *val = '\0';
1349 ++val;
1350
1351 POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1352
1353 if ( !strcmp(key, "join") )
1354 {
1355 if ( !strcmp(val, "round") )
1356 {
1357 joinStyle = JOIN_ROUND;
1358 }
1359 else if ( !(strcmp(val, "mitre") && strcmp(val, "miter")) )
1360 {
1361 joinStyle = JOIN_MITRE;
1362 }
1363 else if ( ! strcmp(val, "bevel") )
1364 {
1365 joinStyle = JOIN_BEVEL;
1366 }
1367 else
1368 {
1369 lwpgerror(
1370 "Invalid buffer end cap style: %s (accept: 'round', 'mitre', 'miter' or 'bevel')",
1371 val);
1372 break;
1373 }
1374 }
1375 else if ( !strcmp(key, "mitre_limit") ||
1376 !strcmp(key, "miter_limit") )
1377 {
1378 /* mitreLimit is a float */
1379 mitreLimit = atof(val);
1380 }
1381 else if ( !strcmp(key, "quad_segs") )
1382 {
1383 /* quadrant segments is an int */
1384 quadsegs = atoi(val);
1385 }
1386 else
1387 {
1388 lwpgerror(
1389 "Invalid buffer parameter: %s (accept: 'join', 'mitre_limit', 'miter_limit and 'quad_segs')",
1390 key);
1391 break;
1392 }
1393 }
1394 POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g", joinStyle, mitreLimit);
1395 pfree(paramstr); /* alloc'ed in text_to_cstring */
1396 }
1397
1398 lwgeom_result = lwgeom_offsetcurve(lwgeom_input, size, quadsegs, joinStyle, mitreLimit);
1399
1400 if (!lwgeom_result)
1401 lwpgerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");
1402
1403 gser_result = geometry_serialize(lwgeom_result);
1404 lwgeom_free(lwgeom_input);
1405 lwgeom_free(lwgeom_result);
1406 PG_RETURN_POINTER(gser_result);
1407 }
1408
1409 PG_FUNCTION_INFO_V1(ST_Intersection);
ST_Intersection(PG_FUNCTION_ARGS)1410 Datum ST_Intersection(PG_FUNCTION_ARGS)
1411 {
1412 GSERIALIZED *geom1;
1413 GSERIALIZED *geom2;
1414 GSERIALIZED *result;
1415 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
1416 double prec = -1;
1417
1418 geom1 = PG_GETARG_GSERIALIZED_P(0);
1419 geom2 = PG_GETARG_GSERIALIZED_P(1);
1420 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
1421 prec = PG_GETARG_FLOAT8(2);
1422
1423 lwgeom1 = lwgeom_from_gserialized(geom1);
1424 lwgeom2 = lwgeom_from_gserialized(geom2);
1425
1426 lwresult = lwgeom_intersection_prec(lwgeom1, lwgeom2, prec);
1427 result = geometry_serialize(lwresult);
1428
1429 lwgeom_free(lwgeom1);
1430 lwgeom_free(lwgeom2);
1431 lwgeom_free(lwresult);
1432
1433 PG_FREE_IF_COPY(geom1, 0);
1434 PG_FREE_IF_COPY(geom2, 1);
1435
1436 PG_RETURN_POINTER(result);
1437 }
1438
1439 PG_FUNCTION_INFO_V1(ST_Difference);
ST_Difference(PG_FUNCTION_ARGS)1440 Datum ST_Difference(PG_FUNCTION_ARGS)
1441 {
1442 GSERIALIZED *geom1;
1443 GSERIALIZED *geom2;
1444 GSERIALIZED *result;
1445 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
1446 double prec = -1;
1447
1448 geom1 = PG_GETARG_GSERIALIZED_P(0);
1449 geom2 = PG_GETARG_GSERIALIZED_P(1);
1450 if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
1451 prec = PG_GETARG_FLOAT8(2);
1452
1453 lwgeom1 = lwgeom_from_gserialized(geom1);
1454 lwgeom2 = lwgeom_from_gserialized(geom2);
1455
1456 lwresult = lwgeom_difference_prec(lwgeom1, lwgeom2, prec);
1457 result = geometry_serialize(lwresult);
1458
1459 lwgeom_free(lwgeom1);
1460 lwgeom_free(lwgeom2);
1461 lwgeom_free(lwresult);
1462
1463 PG_FREE_IF_COPY(geom1, 0);
1464 PG_FREE_IF_COPY(geom2, 1);
1465
1466 PG_RETURN_POINTER(result);
1467 }
1468
1469 /**
1470 @example pointonsurface - {@link #pointonsurface} SELECT ST_PointOnSurface('POLYGON((0 0, 10 0, 10 10, 0 10, 0
1471 0))');
1472 */
1473 PG_FUNCTION_INFO_V1(pointonsurface);
pointonsurface(PG_FUNCTION_ARGS)1474 Datum pointonsurface(PG_FUNCTION_ARGS)
1475 {
1476 GSERIALIZED *geom, *result;
1477 LWGEOM *lwgeom, *lwresult;
1478
1479 geom = PG_GETARG_GSERIALIZED_P(0);
1480
1481 lwgeom = lwgeom_from_gserialized(geom);
1482 lwresult = lwgeom_pointonsurface(lwgeom);
1483 lwgeom_free(lwgeom);
1484 PG_FREE_IF_COPY(geom, 0);
1485
1486 if (!lwresult) PG_RETURN_NULL();
1487
1488 result = geometry_serialize(lwresult);
1489 lwgeom_free(lwresult);
1490 PG_RETURN_POINTER(result);
1491 }
1492
1493 PG_FUNCTION_INFO_V1(centroid);
centroid(PG_FUNCTION_ARGS)1494 Datum centroid(PG_FUNCTION_ARGS)
1495 {
1496 GSERIALIZED *geom, *result;
1497 LWGEOM *lwgeom, *lwresult;
1498
1499 geom = PG_GETARG_GSERIALIZED_P(0);
1500
1501 lwgeom = lwgeom_from_gserialized(geom);
1502 lwresult = lwgeom_centroid(lwgeom);
1503 lwgeom_free(lwgeom);
1504 PG_FREE_IF_COPY(geom, 0);
1505
1506 if (!lwresult) PG_RETURN_NULL();
1507
1508 result = geometry_serialize(lwresult);
1509 lwgeom_free(lwresult);
1510 PG_RETURN_POINTER(result);
1511 }
1512
1513 Datum ST_ReducePrecision(PG_FUNCTION_ARGS);
1514 PG_FUNCTION_INFO_V1(ST_ReducePrecision);
ST_ReducePrecision(PG_FUNCTION_ARGS)1515 Datum ST_ReducePrecision(PG_FUNCTION_ARGS)
1516 {
1517 GSERIALIZED *geom, *result;
1518 LWGEOM *lwgeom, *lwresult;
1519 double gridSize = PG_GETARG_FLOAT8(1);
1520 geom = PG_GETARG_GSERIALIZED_P(0);
1521
1522 lwgeom = lwgeom_from_gserialized(geom);
1523 lwresult = lwgeom_reduceprecision(lwgeom, gridSize);
1524 lwgeom_free(lwgeom);
1525 PG_FREE_IF_COPY(geom, 0);
1526
1527 if (!lwresult) PG_RETURN_NULL();
1528
1529 result = geometry_serialize(lwresult);
1530 lwgeom_free(lwresult);
1531 PG_RETURN_POINTER(result);
1532 }
1533
1534 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS);
1535 PG_FUNCTION_INFO_V1(ST_ClipByBox2d);
ST_ClipByBox2d(PG_FUNCTION_ARGS)1536 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS)
1537 {
1538 static const uint32_t geom_idx = 0;
1539 static const uint32_t box2d_idx = 1;
1540 GSERIALIZED *result;
1541 LWGEOM *lwgeom1, *lwresult ;
1542 GBOX bbox1 = {0};
1543 GBOX *bbox2;
1544 uint8_t type;
1545 int32_t srid;
1546 lwflags_t flags;
1547
1548 if (!gserialized_datum_get_internals_p(PG_GETARG_DATUM(geom_idx), &bbox1, &flags, &type, &srid))
1549 {
1550 /* empty clips to empty, no matter rect */
1551 PG_RETURN_DATUM(PG_GETARG_DATUM(geom_idx));
1552 }
1553
1554 /* WARNING: this is really a BOX2DF, use only xmin and ymin fields */
1555 bbox2 = (GBOX *)PG_GETARG_POINTER(box2d_idx);
1556 bbox2->flags = 0;
1557
1558 /* if bbox1 is covered by bbox2, return lwgeom1 */
1559 if (gbox_contains_2d(bbox2, &bbox1))
1560 {
1561 PG_RETURN_DATUM(PG_GETARG_DATUM(geom_idx));
1562 }
1563
1564 /* If bbox1 outside of bbox2, return empty */
1565 if (!gbox_overlaps_2d(&bbox1, bbox2))
1566 {
1567 /* Get type and srid from datum */
1568 lwresult = lwgeom_construct_empty(type, srid, 0, 0);
1569 result = geometry_serialize(lwresult) ;
1570 lwgeom_free(lwresult) ;
1571 PG_RETURN_POINTER(result);
1572 }
1573
1574 lwgeom1 = lwgeom_from_gserialized(PG_GETARG_GSERIALIZED_P(geom_idx));
1575 lwresult = lwgeom_clip_by_rect(lwgeom1, bbox2->xmin, bbox2->ymin,
1576 bbox2->xmax, bbox2->ymax);
1577
1578 lwgeom_free(lwgeom1);
1579
1580 if (!lwresult)
1581 PG_RETURN_NULL();
1582
1583 result = geometry_serialize(lwresult) ;
1584 PG_RETURN_POINTER(result);
1585 }
1586
1587
1588 /*---------------------------------------------*/
1589
1590 PG_FUNCTION_INFO_V1(isvalid);
isvalid(PG_FUNCTION_ARGS)1591 Datum isvalid(PG_FUNCTION_ARGS)
1592 {
1593 GSERIALIZED *geom1;
1594 LWGEOM *lwgeom;
1595 char result;
1596 GEOSGeom g1;
1597
1598 geom1 = PG_GETARG_GSERIALIZED_P(0);
1599
1600 /* Empty.IsValid() == TRUE */
1601 if ( gserialized_is_empty(geom1) )
1602 PG_RETURN_BOOL(true);
1603
1604 initGEOS(lwpgnotice, lwgeom_geos_error);
1605
1606 lwgeom = lwgeom_from_gserialized(geom1);
1607 if ( ! lwgeom )
1608 {
1609 lwpgerror("unable to deserialize input");
1610 }
1611 g1 = LWGEOM2GEOS(lwgeom, 0);
1612 lwgeom_free(lwgeom);
1613
1614 if ( ! g1 )
1615 {
1616 /* should we drop the following
1617 * notice now that we have ST_isValidReason ?
1618 */
1619 lwpgnotice("%s", lwgeom_geos_errmsg);
1620 PG_RETURN_BOOL(false);
1621 }
1622
1623 result = GEOSisValid(g1);
1624 GEOSGeom_destroy(g1);
1625
1626 if (result == 2)
1627 {
1628 elog(ERROR,"GEOS isvalid() threw an error!");
1629 PG_RETURN_NULL(); /*never get here */
1630 }
1631
1632 PG_FREE_IF_COPY(geom1, 0);
1633 PG_RETURN_BOOL(result);
1634 }
1635
1636 PG_FUNCTION_INFO_V1(isvalidreason);
isvalidreason(PG_FUNCTION_ARGS)1637 Datum isvalidreason(PG_FUNCTION_ARGS)
1638 {
1639 GSERIALIZED *geom = NULL;
1640 char *reason_str = NULL;
1641 text *result = NULL;
1642 const GEOSGeometry *g1 = NULL;
1643
1644 geom = PG_GETARG_GSERIALIZED_P(0);
1645
1646 initGEOS(lwpgnotice, lwgeom_geos_error);
1647
1648 g1 = POSTGIS2GEOS(geom);
1649 if ( g1 )
1650 {
1651 reason_str = GEOSisValidReason(g1);
1652 GEOSGeom_destroy((GEOSGeometry *)g1);
1653 if (!reason_str) HANDLE_GEOS_ERROR("GEOSisValidReason");
1654 result = cstring_to_text(reason_str);
1655 GEOSFree(reason_str);
1656 }
1657 else
1658 {
1659 result = cstring_to_text(lwgeom_geos_errmsg);
1660 }
1661
1662 PG_FREE_IF_COPY(geom, 0);
1663 PG_RETURN_POINTER(result);
1664 }
1665
1666 PG_FUNCTION_INFO_V1(isvaliddetail);
isvaliddetail(PG_FUNCTION_ARGS)1667 Datum isvaliddetail(PG_FUNCTION_ARGS)
1668 {
1669 GSERIALIZED *geom = NULL;
1670 const GEOSGeometry *g1 = NULL;
1671 char *values[3]; /* valid bool, reason text, location geometry */
1672 char *geos_reason = NULL;
1673 char *reason = NULL;
1674 GEOSGeometry *geos_location = NULL;
1675 LWGEOM *location = NULL;
1676 char valid = 0;
1677 HeapTupleHeader result;
1678 TupleDesc tupdesc;
1679 HeapTuple tuple;
1680 AttInMetadata *attinmeta;
1681 int flags = 0;
1682
1683 /*
1684 * Build a tuple description for a
1685 * valid_detail tuple
1686 */
1687 get_call_result_type(fcinfo, 0, &tupdesc);
1688 BlessTupleDesc(tupdesc);
1689
1690 /*
1691 * generate attribute metadata needed later to produce
1692 * tuples from raw C strings
1693 */
1694 attinmeta = TupleDescGetAttInMetadata(tupdesc);
1695
1696 geom = PG_GETARG_GSERIALIZED_P(0);
1697 flags = PG_GETARG_INT32(1);
1698
1699 initGEOS(lwpgnotice, lwgeom_geos_error);
1700
1701 g1 = POSTGIS2GEOS(geom);
1702
1703 if ( g1 )
1704 {
1705 valid = GEOSisValidDetail(g1, flags, &geos_reason, &geos_location);
1706 GEOSGeom_destroy((GEOSGeometry *)g1);
1707 if ( geos_reason )
1708 {
1709 reason = pstrdup(geos_reason);
1710 GEOSFree(geos_reason);
1711 }
1712 if ( geos_location )
1713 {
1714 location = GEOS2LWGEOM(geos_location, GEOSHasZ(geos_location));
1715 GEOSGeom_destroy(geos_location);
1716 }
1717
1718 if (valid == 2)
1719 {
1720 /* NOTE: should only happen on OOM or similar */
1721 lwpgerror("GEOS isvaliddetail() threw an exception!");
1722 PG_RETURN_NULL(); /* never gets here */
1723 }
1724 }
1725 else
1726 {
1727 /* TODO: check lwgeom_geos_errmsg for validity error */
1728 reason = pstrdup(lwgeom_geos_errmsg);
1729 }
1730
1731 /* the boolean validity */
1732 values[0] = valid ? "t" : "f";
1733
1734 /* the reason */
1735 values[1] = reason;
1736
1737 /* the location */
1738 values[2] = location ? lwgeom_to_hexwkb_buffer(location, WKB_EXTENDED) : 0;
1739
1740 tuple = BuildTupleFromCStrings(attinmeta, values);
1741 result = (HeapTupleHeader) palloc(tuple->t_len);
1742 memcpy(result, tuple->t_data, tuple->t_len);
1743 heap_freetuple(tuple);
1744
1745 PG_RETURN_HEAPTUPLEHEADER(result);
1746 }
1747
1748 /**
1749 * overlaps(GSERIALIZED g1,GSERIALIZED g2)
1750 * @param g1
1751 * @param g2
1752 * @return if GEOS::g1->overlaps(g2) returns true
1753 * @throw an error (elog(ERROR,...)) if GEOS throws an error
1754 */
1755 PG_FUNCTION_INFO_V1(overlaps);
overlaps(PG_FUNCTION_ARGS)1756 Datum overlaps(PG_FUNCTION_ARGS)
1757 {
1758 GSERIALIZED *geom1;
1759 GSERIALIZED *geom2;
1760 GEOSGeometry *g1, *g2;
1761 char result;
1762 GBOX box1, box2;
1763
1764 geom1 = PG_GETARG_GSERIALIZED_P(0);
1765 geom2 = PG_GETARG_GSERIALIZED_P(1);
1766 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1767
1768 /* A.Overlaps(Empty) == FALSE */
1769 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1770 PG_RETURN_BOOL(false);
1771
1772 /*
1773 * short-circuit 1: if geom2 bounding box does not overlap
1774 * geom1 bounding box we can return FALSE.
1775 */
1776 if ( gserialized_get_gbox_p(geom1, &box1) &&
1777 gserialized_get_gbox_p(geom2, &box2) )
1778 {
1779 if ( ! gbox_overlaps_2d(&box1, &box2) )
1780 {
1781 PG_RETURN_BOOL(false);
1782 }
1783 }
1784
1785 initGEOS(lwpgnotice, lwgeom_geos_error);
1786
1787 g1 = POSTGIS2GEOS(geom1);
1788 if (!g1)
1789 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1790
1791 g2 = POSTGIS2GEOS(geom2);
1792
1793 if (!g2)
1794 {
1795 GEOSGeom_destroy(g1);
1796 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1797 }
1798
1799 result = GEOSOverlaps(g1,g2);
1800
1801 GEOSGeom_destroy(g1);
1802 GEOSGeom_destroy(g2);
1803 if (result == 2) HANDLE_GEOS_ERROR("GEOSOverlaps");
1804
1805 PG_FREE_IF_COPY(geom1, 0);
1806 PG_FREE_IF_COPY(geom2, 1);
1807
1808 PG_RETURN_BOOL(result);
1809 }
1810
1811
1812 PG_FUNCTION_INFO_V1(contains);
contains(PG_FUNCTION_ARGS)1813 Datum contains(PG_FUNCTION_ARGS)
1814 {
1815 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
1816 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
1817 const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
1818 const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
1819 int result;
1820 GEOSGeometry *g1, *g2;
1821 GBOX box1, box2;
1822 PrepGeomCache *prep_cache;
1823 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1824
1825 /* A.Contains(Empty) == FALSE */
1826 if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
1827 PG_RETURN_BOOL(false);
1828
1829 POSTGIS_DEBUG(3, "contains called.");
1830
1831 /*
1832 ** short-circuit 1: if geom2 bounding box is not completely inside
1833 ** geom1 bounding box we can return FALSE.
1834 */
1835 if (gserialized_get_gbox_p(geom1, &box1) &&
1836 gserialized_get_gbox_p(geom2, &box2))
1837 {
1838 if (!gbox_contains_2d(&box1, &box2))
1839 PG_RETURN_BOOL(false);
1840 }
1841
1842 /*
1843 ** short-circuit 2: if geom2 is a point and geom1 is a polygon
1844 ** call the point-in-polygon function.
1845 */
1846 if (is_poly(geom1) && is_point(geom2))
1847 {
1848 SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
1849 SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
1850 const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
1851 const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
1852 RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
1853 int retval;
1854
1855 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1856 if (gserialized_get_type(gpoint) == POINTTYPE)
1857 {
1858 LWGEOM* point = lwgeom_from_gserialized(gpoint);
1859 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1860 lwgeom_free(point);
1861
1862 retval = (pip_result == 1); /* completely inside */
1863 }
1864 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
1865 {
1866 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
1867 uint32_t i;
1868 int found_completely_inside = LW_FALSE;
1869
1870 retval = LW_TRUE;
1871 for (i = 0; i < mpoint->ngeoms; i++)
1872 {
1873 /* We need to find at least one point that's completely inside the
1874 * polygons (pip_result == 1). As long as we have one point that's
1875 * completely inside, we can have as many as we want on the boundary
1876 * itself. (pip_result == 0)
1877 */
1878 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
1879 if (pip_result == 1)
1880 found_completely_inside = LW_TRUE;
1881
1882 if (pip_result == -1) /* completely outside */
1883 {
1884 retval = LW_FALSE;
1885 break;
1886 }
1887 }
1888
1889 retval = retval && found_completely_inside;
1890 lwmpoint_free(mpoint);
1891 }
1892 else
1893 {
1894 /* Never get here */
1895 elog(ERROR,"Type isn't point or multipoint!");
1896 PG_RETURN_BOOL(false);
1897 }
1898
1899 return retval > 0;
1900 }
1901 else
1902 {
1903 POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
1904 }
1905
1906 initGEOS(lwpgnotice, lwgeom_geos_error);
1907
1908 prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, NULL);
1909
1910 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1911 {
1912 g1 = POSTGIS2GEOS(geom2);
1913 if (!g1)
1914 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
1915
1916 POSTGIS_DEBUG(4, "containsPrepared: cache is live, running preparedcontains");
1917 result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
1918 GEOSGeom_destroy(g1);
1919 }
1920 else
1921 {
1922 g1 = POSTGIS2GEOS(geom1);
1923 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1924 g2 = POSTGIS2GEOS(geom2);
1925 if (!g2)
1926 {
1927 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1928 GEOSGeom_destroy(g1);
1929 }
1930 POSTGIS_DEBUG(4, "containsPrepared: cache is not ready, running standard contains");
1931 result = GEOSContains( g1, g2);
1932 GEOSGeom_destroy(g1);
1933 GEOSGeom_destroy(g2);
1934 }
1935
1936 if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1937
1938 PG_RETURN_BOOL(result > 0);
1939 }
1940
1941
1942 PG_FUNCTION_INFO_V1(within);
within(PG_FUNCTION_ARGS)1943 Datum within(PG_FUNCTION_ARGS)
1944 {
1945 PG_RETURN_DATUM(CallerFInfoFunctionCall2(contains, fcinfo->flinfo, InvalidOid,
1946 PG_GETARG_DATUM(1), PG_GETARG_DATUM(0)));
1947 }
1948
1949
1950
1951 PG_FUNCTION_INFO_V1(containsproperly);
containsproperly(PG_FUNCTION_ARGS)1952 Datum containsproperly(PG_FUNCTION_ARGS)
1953 {
1954 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
1955 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
1956 const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
1957 const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
1958 char result;
1959 GBOX box1, box2;
1960 PrepGeomCache * prep_cache;
1961
1962 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1963
1964 /* A.ContainsProperly(Empty) == FALSE */
1965 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1966 PG_RETURN_BOOL(false);
1967
1968 /*
1969 * short-circuit: if geom2 bounding box is not completely inside
1970 * geom1 bounding box we can return FALSE.
1971 */
1972 if ( gserialized_get_gbox_p(geom1, &box1) &&
1973 gserialized_get_gbox_p(geom2, &box2) )
1974 {
1975 if ( ! gbox_contains_2d(&box1, &box2) )
1976 PG_RETURN_BOOL(false);
1977 }
1978
1979 initGEOS(lwpgnotice, lwgeom_geos_error);
1980
1981 prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
1982
1983 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1984 {
1985 GEOSGeometry *g = POSTGIS2GEOS(geom2);
1986 if (!g) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1987 result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
1988 GEOSGeom_destroy(g);
1989 }
1990 else
1991 {
1992 GEOSGeometry *g2;
1993 GEOSGeometry *g1;
1994
1995 g1 = POSTGIS2GEOS(geom1);
1996 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1997 g2 = POSTGIS2GEOS(geom2);
1998 if (!g2)
1999 {
2000 GEOSGeom_destroy(g1);
2001 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2002 }
2003 result = GEOSRelatePattern( g1, g2, "T**FF*FF*" );
2004 GEOSGeom_destroy(g1);
2005 GEOSGeom_destroy(g2);
2006 }
2007
2008 if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
2009
2010 PG_RETURN_BOOL(result);
2011 }
2012
2013 /*
2014 * Described at
2015 * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2016 */
2017 PG_FUNCTION_INFO_V1(covers);
covers(PG_FUNCTION_ARGS)2018 Datum covers(PG_FUNCTION_ARGS)
2019 {
2020 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2021 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2022 const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2023 const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2024 int result;
2025 GBOX box1, box2;
2026 PrepGeomCache *prep_cache;
2027
2028
2029 /* A.Covers(Empty) == FALSE */
2030 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2031 PG_RETURN_BOOL(false);
2032
2033 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2034
2035 /*
2036 * short-circuit 1: if geom2 bounding box is not completely inside
2037 * geom1 bounding box we can return FALSE.
2038 */
2039 if ( gserialized_get_gbox_p(geom1, &box1) &&
2040 gserialized_get_gbox_p(geom2, &box2) )
2041 {
2042 if ( ! gbox_contains_2d(&box1, &box2) )
2043 {
2044 PG_RETURN_BOOL(false);
2045 }
2046 }
2047 /*
2048 * short-circuit 2: if geom2 is a point and geom1 is a polygon
2049 * call the point-in-polygon function.
2050 */
2051 if (is_poly(geom1) && is_point(geom2))
2052 {
2053 SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2054 SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2055 const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2056 const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2057 RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2058 int retval;
2059
2060 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2061 if (gserialized_get_type(gpoint) == POINTTYPE)
2062 {
2063 LWGEOM* point = lwgeom_from_gserialized(gpoint);
2064 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2065 lwgeom_free(point);
2066
2067 retval = (pip_result != -1); /* not outside */
2068 }
2069 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2070 {
2071 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2072 uint32_t i;
2073
2074 retval = LW_TRUE;
2075 for (i = 0; i < mpoint->ngeoms; i++)
2076 {
2077 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2078 if (pip_result == -1)
2079 {
2080 retval = LW_FALSE;
2081 break;
2082 }
2083 }
2084
2085 lwmpoint_free(mpoint);
2086 }
2087 else
2088 {
2089 /* Never get here */
2090 elog(ERROR,"Type isn't point or multipoint!");
2091 PG_RETURN_NULL();
2092 }
2093
2094 PG_RETURN_BOOL(retval);
2095 }
2096 else
2097 {
2098 POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2099 }
2100
2101 initGEOS(lwpgnotice, lwgeom_geos_error);
2102
2103 prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
2104
2105 if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
2106 {
2107 GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
2108 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2109 result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
2110 GEOSGeom_destroy(g1);
2111 }
2112 else
2113 {
2114 GEOSGeometry *g1;
2115 GEOSGeometry *g2;
2116
2117 g1 = POSTGIS2GEOS(geom1);
2118 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2119 g2 = POSTGIS2GEOS(geom2);
2120 if (!g2)
2121 {
2122 GEOSGeom_destroy(g1);
2123 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2124 }
2125 result = GEOSRelatePattern( g1, g2, "******FF*" );
2126 GEOSGeom_destroy(g1);
2127 GEOSGeom_destroy(g2);
2128 }
2129
2130 if (result == 2) HANDLE_GEOS_ERROR("GEOSCovers");
2131
2132 PG_RETURN_BOOL(result);
2133 }
2134
2135
2136 /**
2137 * ST_Within(A, B) => ST_Contains(B, A) so we just delegate this calculation to the
2138 * Contains implementation.
2139 PG_FUNCTION_INFO_V1(within);
2140 Datum within(PG_FUNCTION_ARGS)
2141 */
2142
2143 /*
2144 * Described at:
2145 * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2146 */
2147 PG_FUNCTION_INFO_V1(coveredby);
coveredby(PG_FUNCTION_ARGS)2148 Datum coveredby(PG_FUNCTION_ARGS)
2149 {
2150 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2151 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2152 const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2153 const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2154 GEOSGeometry *g1, *g2;
2155 int result;
2156 GBOX box1, box2;
2157 char *patt = "**F**F***";
2158
2159 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2160
2161 /* A.CoveredBy(Empty) == FALSE */
2162 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2163 PG_RETURN_BOOL(false);
2164
2165 /*
2166 * short-circuit 1: if geom1 bounding box is not completely inside
2167 * geom2 bounding box we can return FALSE.
2168 */
2169 if ( gserialized_get_gbox_p(geom1, &box1) &&
2170 gserialized_get_gbox_p(geom2, &box2) )
2171 {
2172 if ( ! gbox_contains_2d(&box2, &box1) )
2173 {
2174 PG_RETURN_BOOL(false);
2175 }
2176
2177 POSTGIS_DEBUG(3, "bounding box short-circuit missed.");
2178 }
2179 /*
2180 * short-circuit 2: if geom1 is a point and geom2 is a polygon
2181 * call the point-in-polygon function.
2182 */
2183 if (is_point(geom1) && is_poly(geom2))
2184 {
2185 SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2186 SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2187 const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2188 const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2189 RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2190 int retval;
2191
2192 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2193 if (gserialized_get_type(gpoint) == POINTTYPE)
2194 {
2195 LWGEOM* point = lwgeom_from_gserialized(gpoint);
2196 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2197 lwgeom_free(point);
2198
2199 retval = (pip_result != -1); /* not outside */
2200 }
2201 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2202 {
2203 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2204 uint32_t i;
2205
2206 retval = LW_TRUE;
2207 for (i = 0; i < mpoint->ngeoms; i++)
2208 {
2209 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2210 if (pip_result == -1)
2211 {
2212 retval = LW_FALSE;
2213 break;
2214 }
2215 }
2216
2217 lwmpoint_free(mpoint);
2218 }
2219 else
2220 {
2221 /* Never get here */
2222 elog(ERROR,"Type isn't point or multipoint!");
2223 PG_RETURN_NULL();
2224 }
2225
2226 PG_RETURN_BOOL(retval);
2227 }
2228 else
2229 {
2230 POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2231 }
2232
2233 initGEOS(lwpgnotice, lwgeom_geos_error);
2234
2235 g1 = POSTGIS2GEOS(geom1);
2236
2237 if (!g1)
2238 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2239
2240 g2 = POSTGIS2GEOS(geom2);
2241
2242 if (!g2)
2243 {
2244 GEOSGeom_destroy(g1);
2245 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2246 }
2247
2248 result = GEOSRelatePattern(g1,g2,patt);
2249
2250 GEOSGeom_destroy(g1);
2251 GEOSGeom_destroy(g2);
2252
2253 if (result == 2) HANDLE_GEOS_ERROR("GEOSCoveredBy");
2254
2255 PG_RETURN_BOOL(result);
2256 }
2257
2258 PG_FUNCTION_INFO_V1(crosses);
crosses(PG_FUNCTION_ARGS)2259 Datum crosses(PG_FUNCTION_ARGS)
2260 {
2261 GSERIALIZED *geom1;
2262 GSERIALIZED *geom2;
2263 GEOSGeometry *g1, *g2;
2264 int result;
2265 GBOX box1, box2;
2266
2267 geom1 = PG_GETARG_GSERIALIZED_P(0);
2268 geom2 = PG_GETARG_GSERIALIZED_P(1);
2269 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2270
2271 /* A.Crosses(Empty) == FALSE */
2272 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2273 PG_RETURN_BOOL(false);
2274
2275 /*
2276 * short-circuit 1: if geom2 bounding box does not overlap
2277 * geom1 bounding box we can return FALSE.
2278 */
2279 if ( gserialized_get_gbox_p(geom1, &box1) &&
2280 gserialized_get_gbox_p(geom2, &box2) )
2281 {
2282 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2283 {
2284 PG_RETURN_BOOL(false);
2285 }
2286 }
2287
2288 initGEOS(lwpgnotice, lwgeom_geos_error);
2289
2290 g1 = POSTGIS2GEOS(geom1);
2291 if (!g1)
2292 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2293
2294 g2 = POSTGIS2GEOS(geom2);
2295 if (!g2)
2296 {
2297 GEOSGeom_destroy(g1);
2298 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2299 }
2300
2301 result = GEOSCrosses(g1,g2);
2302
2303 GEOSGeom_destroy(g1);
2304 GEOSGeom_destroy(g2);
2305
2306 if (result == 2) HANDLE_GEOS_ERROR("GEOSCrosses");
2307
2308 PG_FREE_IF_COPY(geom1, 0);
2309 PG_FREE_IF_COPY(geom2, 1);
2310
2311 PG_RETURN_BOOL(result);
2312 }
2313
2314 PG_FUNCTION_INFO_V1(ST_Intersects);
ST_Intersects(PG_FUNCTION_ARGS)2315 Datum ST_Intersects(PG_FUNCTION_ARGS)
2316 {
2317 SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2318 SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2319 const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2320 const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2321 int result;
2322 GBOX box1, box2;
2323 PrepGeomCache *prep_cache;
2324
2325 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2326
2327 /* A.Intersects(Empty) == FALSE */
2328 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2329 PG_RETURN_BOOL(false);
2330
2331 /*
2332 * short-circuit 1: if geom2 bounding box does not overlap
2333 * geom1 bounding box we can return FALSE.
2334 */
2335 if ( gserialized_get_gbox_p(geom1, &box1) &&
2336 gserialized_get_gbox_p(geom2, &box2) )
2337 {
2338 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2339 PG_RETURN_BOOL(false);
2340 }
2341
2342 /*
2343 * short-circuit 2: if the geoms are a point and a polygon,
2344 * call the point_outside_polygon function.
2345 */
2346 if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
2347 {
2348 SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2349 SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2350 const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2351 const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2352 RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2353 int retval;
2354
2355 POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2356 if (gserialized_get_type(gpoint) == POINTTYPE)
2357 {
2358 LWGEOM* point = lwgeom_from_gserialized(gpoint);
2359 int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2360 lwgeom_free(point);
2361
2362 retval = (pip_result != -1); /* not outside */
2363 }
2364 else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2365 {
2366 LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2367 uint32_t i;
2368
2369 retval = LW_FALSE;
2370 for (i = 0; i < mpoint->ngeoms; i++)
2371 {
2372 int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2373 if (pip_result != -1) /* not outside */
2374 {
2375 retval = LW_TRUE;
2376 break;
2377 }
2378 }
2379
2380 lwmpoint_free(mpoint);
2381 }
2382 else
2383 {
2384 /* Never get here */
2385 elog(ERROR,"Type isn't point or multipoint!");
2386 PG_RETURN_NULL();
2387 }
2388
2389 PG_RETURN_BOOL(retval);
2390 }
2391
2392 initGEOS(lwpgnotice, lwgeom_geos_error);
2393 prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
2394
2395 if ( prep_cache && prep_cache->prepared_geom )
2396 {
2397 if ( prep_cache->gcache.argnum == 1 )
2398 {
2399 GEOSGeometry *g = POSTGIS2GEOS(geom2);
2400 if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2401 result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2402 GEOSGeom_destroy(g);
2403 }
2404 else
2405 {
2406 GEOSGeometry *g = POSTGIS2GEOS(geom1);
2407 if (!g)
2408 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2409 result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2410 GEOSGeom_destroy(g);
2411 }
2412 }
2413 else
2414 {
2415 GEOSGeometry *g1;
2416 GEOSGeometry *g2;
2417 g1 = POSTGIS2GEOS(geom1);
2418 if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2419 g2 = POSTGIS2GEOS(geom2);
2420 if (!g2)
2421 {
2422 GEOSGeom_destroy(g1);
2423 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2424 }
2425 result = GEOSIntersects( g1, g2);
2426 GEOSGeom_destroy(g1);
2427 GEOSGeom_destroy(g2);
2428 }
2429
2430 if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
2431
2432 PG_RETURN_BOOL(result);
2433 }
2434
2435
2436 PG_FUNCTION_INFO_V1(touches);
touches(PG_FUNCTION_ARGS)2437 Datum touches(PG_FUNCTION_ARGS)
2438 {
2439 GSERIALIZED *geom1;
2440 GSERIALIZED *geom2;
2441 GEOSGeometry *g1, *g2;
2442 char result;
2443 GBOX box1, box2;
2444
2445 geom1 = PG_GETARG_GSERIALIZED_P(0);
2446 geom2 = PG_GETARG_GSERIALIZED_P(1);
2447 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2448
2449 /* A.Touches(Empty) == FALSE */
2450 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2451 PG_RETURN_BOOL(false);
2452
2453 /*
2454 * short-circuit 1: if geom2 bounding box does not overlap
2455 * geom1 bounding box we can return FALSE.
2456 */
2457 if ( gserialized_get_gbox_p(geom1, &box1) &&
2458 gserialized_get_gbox_p(geom2, &box2) )
2459 {
2460 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2461 {
2462 PG_RETURN_BOOL(false);
2463 }
2464 }
2465
2466 initGEOS(lwpgnotice, lwgeom_geos_error);
2467
2468 g1 = POSTGIS2GEOS(geom1 );
2469 if (!g1)
2470 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2471
2472 g2 = POSTGIS2GEOS(geom2 );
2473 if (!g2)
2474 {
2475 GEOSGeom_destroy(g1);
2476 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2477 }
2478
2479 result = GEOSTouches(g1,g2);
2480
2481 GEOSGeom_destroy(g1);
2482 GEOSGeom_destroy(g2);
2483
2484 if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
2485
2486 PG_FREE_IF_COPY(geom1, 0);
2487 PG_FREE_IF_COPY(geom2, 1);
2488
2489 PG_RETURN_BOOL(result);
2490 }
2491
2492
2493 PG_FUNCTION_INFO_V1(disjoint);
disjoint(PG_FUNCTION_ARGS)2494 Datum disjoint(PG_FUNCTION_ARGS)
2495 {
2496 GSERIALIZED *geom1;
2497 GSERIALIZED *geom2;
2498 GEOSGeometry *g1, *g2;
2499 char result;
2500 GBOX box1, box2;
2501
2502 geom1 = PG_GETARG_GSERIALIZED_P(0);
2503 geom2 = PG_GETARG_GSERIALIZED_P(1);
2504 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2505
2506 /* A.Disjoint(Empty) == TRUE */
2507 if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2508 PG_RETURN_BOOL(true);
2509
2510 /*
2511 * short-circuit 1: if geom2 bounding box does not overlap
2512 * geom1 bounding box we can return TRUE.
2513 */
2514 if ( gserialized_get_gbox_p(geom1, &box1) &&
2515 gserialized_get_gbox_p(geom2, &box2) )
2516 {
2517 if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2518 {
2519 PG_RETURN_BOOL(true);
2520 }
2521 }
2522
2523 initGEOS(lwpgnotice, lwgeom_geos_error);
2524
2525 g1 = POSTGIS2GEOS(geom1);
2526 if (!g1)
2527 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2528
2529 g2 = POSTGIS2GEOS(geom2);
2530 if (!g2)
2531 {
2532 GEOSGeom_destroy(g1);
2533 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2534 }
2535
2536 result = GEOSDisjoint(g1,g2);
2537
2538 GEOSGeom_destroy(g1);
2539 GEOSGeom_destroy(g2);
2540
2541 if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
2542
2543 PG_FREE_IF_COPY(geom1, 0);
2544 PG_FREE_IF_COPY(geom2, 1);
2545
2546 PG_RETURN_BOOL(result);
2547 }
2548
2549
2550 PG_FUNCTION_INFO_V1(relate_pattern);
relate_pattern(PG_FUNCTION_ARGS)2551 Datum relate_pattern(PG_FUNCTION_ARGS)
2552 {
2553 GSERIALIZED *geom1;
2554 GSERIALIZED *geom2;
2555 char *patt;
2556 char result;
2557 GEOSGeometry *g1, *g2;
2558 size_t i;
2559
2560 geom1 = PG_GETARG_GSERIALIZED_P(0);
2561 geom2 = PG_GETARG_GSERIALIZED_P(1);
2562 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2563
2564 /* TODO handle empty */
2565
2566 initGEOS(lwpgnotice, lwgeom_geos_error);
2567
2568 g1 = POSTGIS2GEOS(geom1);
2569 if (!g1)
2570 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2571 g2 = POSTGIS2GEOS(geom2);
2572 if (!g2)
2573 {
2574 GEOSGeom_destroy(g1);
2575 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2576 }
2577
2578 patt = DatumGetCString(DirectFunctionCall1(textout,
2579 PointerGetDatum(PG_GETARG_DATUM(2))));
2580
2581 /*
2582 ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2583 */
2584 for ( i = 0; i < strlen(patt); i++ )
2585 {
2586 if ( patt[i] == 't' ) patt[i] = 'T';
2587 if ( patt[i] == 'f' ) patt[i] = 'F';
2588 }
2589
2590 result = GEOSRelatePattern(g1,g2,patt);
2591 GEOSGeom_destroy(g1);
2592 GEOSGeom_destroy(g2);
2593 pfree(patt);
2594
2595 if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
2596
2597 PG_FREE_IF_COPY(geom1, 0);
2598 PG_FREE_IF_COPY(geom2, 1);
2599
2600 PG_RETURN_BOOL(result);
2601 }
2602
2603
2604
2605 PG_FUNCTION_INFO_V1(relate_full);
relate_full(PG_FUNCTION_ARGS)2606 Datum relate_full(PG_FUNCTION_ARGS)
2607 {
2608 GSERIALIZED *geom1;
2609 GSERIALIZED *geom2;
2610 GEOSGeometry *g1, *g2;
2611 char *relate_str;
2612 text *result;
2613 int bnr = GEOSRELATE_BNR_OGC;
2614
2615 /* TODO handle empty */
2616 geom1 = PG_GETARG_GSERIALIZED_P(0);
2617 geom2 = PG_GETARG_GSERIALIZED_P(1);
2618 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2619
2620 if ( PG_NARGS() > 2 )
2621 bnr = PG_GETARG_INT32(2);
2622
2623 initGEOS(lwpgnotice, lwgeom_geos_error);
2624
2625 g1 = POSTGIS2GEOS(geom1 );
2626 if (!g1)
2627 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2628 g2 = POSTGIS2GEOS(geom2 );
2629 if (!g2)
2630 {
2631 GEOSGeom_destroy(g1);
2632 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2633 }
2634
2635 POSTGIS_DEBUG(3, "constructed geometries ");
2636
2637 POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2638 POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2639
2640 relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2641
2642 GEOSGeom_destroy(g1);
2643 GEOSGeom_destroy(g2);
2644
2645 if (!relate_str) HANDLE_GEOS_ERROR("GEOSRelate");
2646
2647 result = cstring_to_text(relate_str);
2648 GEOSFree(relate_str);
2649
2650 PG_FREE_IF_COPY(geom1, 0);
2651 PG_FREE_IF_COPY(geom2, 1);
2652
2653 PG_RETURN_TEXT_P(result);
2654 }
2655
2656
2657 PG_FUNCTION_INFO_V1(ST_Equals);
ST_Equals(PG_FUNCTION_ARGS)2658 Datum ST_Equals(PG_FUNCTION_ARGS)
2659 {
2660 GSERIALIZED *geom1;
2661 GSERIALIZED *geom2;
2662 GEOSGeometry *g1, *g2;
2663 char result;
2664 GBOX box1, box2;
2665
2666 geom1 = PG_GETARG_GSERIALIZED_P(0);
2667 geom2 = PG_GETARG_GSERIALIZED_P(1);
2668 gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2669
2670 /* Empty == Empty */
2671 if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
2672 PG_RETURN_BOOL(true);
2673
2674 /*
2675 * short-circuit: If geom1 and geom2 do not have the same bounding box
2676 * we can return FALSE.
2677 */
2678 if ( gserialized_get_gbox_p(geom1, &box1) &&
2679 gserialized_get_gbox_p(geom2, &box2) )
2680 {
2681 if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
2682 {
2683 PG_RETURN_BOOL(false);
2684 }
2685 }
2686
2687 /*
2688 * short-circuit: if geom1 and geom2 are binary-equivalent, we can return
2689 * TRUE. This is much faster than doing the comparison using GEOS.
2690 */
2691 if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
2692 PG_RETURN_BOOL(true);
2693 }
2694
2695 initGEOS(lwpgnotice, lwgeom_geos_error);
2696
2697 g1 = POSTGIS2GEOS(geom1);
2698
2699 if (!g1)
2700 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2701
2702 g2 = POSTGIS2GEOS(geom2);
2703
2704 if (!g2)
2705 {
2706 GEOSGeom_destroy(g1);
2707 HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2708 }
2709
2710 result = GEOSEquals(g1,g2);
2711
2712 GEOSGeom_destroy(g1);
2713 GEOSGeom_destroy(g2);
2714
2715 if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
2716
2717 PG_FREE_IF_COPY(geom1, 0);
2718 PG_FREE_IF_COPY(geom2, 1);
2719
2720 PG_RETURN_BOOL(result);
2721 }
2722
2723 PG_FUNCTION_INFO_V1(issimple);
issimple(PG_FUNCTION_ARGS)2724 Datum issimple(PG_FUNCTION_ARGS)
2725 {
2726 GSERIALIZED *geom;
2727 LWGEOM *lwgeom_in;
2728 int result;
2729
2730 POSTGIS_DEBUG(2, "issimple called");
2731
2732 geom = PG_GETARG_GSERIALIZED_P(0);
2733
2734 if ( gserialized_is_empty(geom) )
2735 PG_RETURN_BOOL(true);
2736
2737 lwgeom_in = lwgeom_from_gserialized(geom);
2738 result = lwgeom_is_simple(lwgeom_in);
2739 lwgeom_free(lwgeom_in) ;
2740 PG_FREE_IF_COPY(geom, 0);
2741
2742 if (result == -1) {
2743 PG_RETURN_NULL(); /*never get here */
2744 }
2745
2746 PG_RETURN_BOOL(result);
2747 }
2748
2749 PG_FUNCTION_INFO_V1(isring);
isring(PG_FUNCTION_ARGS)2750 Datum isring(PG_FUNCTION_ARGS)
2751 {
2752 GSERIALIZED *geom;
2753 GEOSGeometry *g1;
2754 int result;
2755
2756 geom = PG_GETARG_GSERIALIZED_P(0);
2757
2758 /* Empty things can't close */
2759 if ( gserialized_is_empty(geom) )
2760 PG_RETURN_BOOL(false);
2761
2762 initGEOS(lwpgnotice, lwgeom_geos_error);
2763
2764 g1 = POSTGIS2GEOS(geom);
2765 if (!g1)
2766 HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2767
2768 if ( GEOSGeomTypeId(g1) != GEOS_LINESTRING )
2769 {
2770 GEOSGeom_destroy(g1);
2771 elog(ERROR, "ST_IsRing() should only be called on a linear feature");
2772 }
2773
2774 result = GEOSisRing(g1);
2775 GEOSGeom_destroy(g1);
2776
2777 if (result == 2) HANDLE_GEOS_ERROR("GEOSisRing");
2778
2779 PG_FREE_IF_COPY(geom, 0);
2780 PG_RETURN_BOOL(result);
2781 }
2782
2783 GSERIALIZED *
GEOS2POSTGIS(GEOSGeom geom,char want3d)2784 GEOS2POSTGIS(GEOSGeom geom, char want3d)
2785 {
2786 LWGEOM *lwgeom;
2787 GSERIALIZED *result;
2788
2789 lwgeom = GEOS2LWGEOM(geom, want3d);
2790 if ( ! lwgeom )
2791 {
2792 lwpgerror("%s: GEOS2LWGEOM returned NULL", __func__);
2793 return NULL;
2794 }
2795
2796 POSTGIS_DEBUGF(4, "%s: GEOS2LWGEOM returned a %s", __func__, lwgeom_summary(lwgeom, 0));
2797
2798 if (lwgeom_needs_bbox(lwgeom)) lwgeom_add_bbox(lwgeom);
2799
2800 result = geometry_serialize(lwgeom);
2801 lwgeom_free(lwgeom);
2802
2803 return result;
2804 }
2805
2806 /*-----=POSTGIS2GEOS= */
2807
2808 GEOSGeometry *
POSTGIS2GEOS(const GSERIALIZED * pglwgeom)2809 POSTGIS2GEOS(const GSERIALIZED *pglwgeom)
2810 {
2811 GEOSGeometry *ret;
2812 LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
2813 if ( ! lwgeom )
2814 {
2815 lwpgerror("POSTGIS2GEOS: unable to deserialize input");
2816 return NULL;
2817 }
2818 ret = LWGEOM2GEOS(lwgeom, 0);
2819 lwgeom_free(lwgeom);
2820
2821 return ret;
2822 }
2823
array_nelems_not_null(ArrayType * array)2824 uint32_t array_nelems_not_null(ArrayType* array) {
2825 ArrayIterator iterator;
2826 Datum value;
2827 bool isnull;
2828 uint32_t nelems_not_null = 0;
2829 iterator = array_create_iterator(array, 0, NULL);
2830 while(array_iterate(iterator, &value, &isnull) )
2831 if (!isnull)
2832 nelems_not_null++;
2833
2834 array_free_iterator(iterator);
2835
2836 return nelems_not_null;
2837 }
2838
2839 /* ARRAY2LWGEOM: Converts the non-null elements of a Postgres array into a LWGEOM* array */
ARRAY2LWGEOM(ArrayType * array,uint32_t nelems,int * is3d,int * srid)2840 LWGEOM** ARRAY2LWGEOM(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2841 {
2842 ArrayIterator iterator;
2843 Datum value;
2844 bool isnull;
2845 bool gotsrid = false;
2846 uint32_t i = 0;
2847
2848 LWGEOM** lw_geoms = palloc(nelems * sizeof(LWGEOM*));
2849
2850 iterator = array_create_iterator(array, 0, NULL);
2851
2852 while (array_iterate(iterator, &value, &isnull))
2853 {
2854 GSERIALIZED *geom = (GSERIALIZED *)DatumGetPointer(value);
2855
2856 if (isnull)
2857 continue;
2858
2859 *is3d = *is3d || gserialized_has_z(geom);
2860
2861 lw_geoms[i] = lwgeom_from_gserialized(geom);
2862 if (!lw_geoms[i]) /* error in creation */
2863 {
2864 lwpgerror("Geometry deserializing geometry");
2865 return NULL;
2866 }
2867 if (!gotsrid)
2868 {
2869 gotsrid = true;
2870 *srid = gserialized_get_srid(geom);
2871 }
2872 else
2873 gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2874
2875 i++;
2876 }
2877
2878 return lw_geoms;
2879 }
2880
2881 /* ARRAY2GEOS: Converts the non-null elements of a Postgres array into a GEOSGeometry* array */
ARRAY2GEOS(ArrayType * array,uint32_t nelems,int * is3d,int * srid)2882 GEOSGeometry** ARRAY2GEOS(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2883 {
2884 ArrayIterator iterator;
2885 Datum value;
2886 bool isnull;
2887 bool gotsrid = false;
2888 uint32_t i = 0;
2889
2890 GEOSGeometry** geos_geoms = palloc(nelems * sizeof(GEOSGeometry*));
2891
2892 iterator = array_create_iterator(array, 0, NULL);
2893
2894 while(array_iterate(iterator, &value, &isnull))
2895 {
2896 GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2897
2898 if (isnull)
2899 continue;
2900
2901 *is3d = *is3d || gserialized_has_z(geom);
2902
2903 geos_geoms[i] = POSTGIS2GEOS(geom);
2904 if (!geos_geoms[i])
2905 {
2906 uint32_t j;
2907 lwpgerror("Geometry could not be converted to GEOS");
2908
2909 for (j = 0; j < i; j++) {
2910 GEOSGeom_destroy(geos_geoms[j]);
2911 }
2912 return NULL;
2913 }
2914
2915 if (!gotsrid)
2916 {
2917 *srid = gserialized_get_srid(geom);
2918 gotsrid = true;
2919 }
2920 else if (*srid != gserialized_get_srid(geom))
2921 {
2922 uint32_t j;
2923 for (j = 0; j <= i; j++) {
2924 GEOSGeom_destroy(geos_geoms[j]);
2925 }
2926 gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2927 return NULL;
2928 }
2929
2930 i++;
2931 }
2932
2933 array_free_iterator(iterator);
2934 return geos_geoms;
2935 }
2936
2937 PG_FUNCTION_INFO_V1(GEOSnoop);
GEOSnoop(PG_FUNCTION_ARGS)2938 Datum GEOSnoop(PG_FUNCTION_ARGS)
2939 {
2940 GSERIALIZED *geom;
2941 GEOSGeometry *geosgeom;
2942 GSERIALIZED *lwgeom_result;
2943
2944 initGEOS(lwpgnotice, lwgeom_geos_error);
2945
2946 geom = PG_GETARG_GSERIALIZED_P(0);
2947 geosgeom = POSTGIS2GEOS(geom);
2948 if ( ! geosgeom ) PG_RETURN_NULL();
2949
2950 lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
2951 GEOSGeom_destroy(geosgeom);
2952
2953 PG_FREE_IF_COPY(geom, 0);
2954
2955 PG_RETURN_POINTER(lwgeom_result);
2956 }
2957
2958 PG_FUNCTION_INFO_V1(polygonize_garray);
polygonize_garray(PG_FUNCTION_ARGS)2959 Datum polygonize_garray(PG_FUNCTION_ARGS)
2960 {
2961 ArrayType *array;
2962 int is3d = 0;
2963 uint32 nelems, i;
2964 GSERIALIZED *result;
2965 GEOSGeometry *geos_result;
2966 const GEOSGeometry **vgeoms;
2967 int32_t srid = SRID_UNKNOWN;
2968 #if POSTGIS_DEBUG_LEVEL >= 3
2969 static int call=1;
2970 #endif
2971
2972 #if POSTGIS_DEBUG_LEVEL >= 3
2973 call++;
2974 #endif
2975
2976 if (PG_ARGISNULL(0))
2977 PG_RETURN_NULL();
2978
2979 array = PG_GETARG_ARRAYTYPE_P(0);
2980 nelems = array_nelems_not_null(array);
2981
2982 if (nelems == 0)
2983 PG_RETURN_NULL();
2984
2985 POSTGIS_DEBUGF(3, "polygonize_garray: number of non-null elements: %d", nelems);
2986
2987 /* Ok, we really need geos now ;) */
2988 initGEOS(lwpgnotice, lwgeom_geos_error);
2989
2990 vgeoms = (const GEOSGeometry**) ARRAY2GEOS(array, nelems, &is3d, &srid);
2991
2992 POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
2993
2994 geos_result = GEOSPolygonize(vgeoms, nelems);
2995
2996 POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
2997
2998 for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
2999 pfree(vgeoms);
3000
3001 if ( ! geos_result ) PG_RETURN_NULL();
3002
3003 GEOSSetSRID(geos_result, srid);
3004 result = GEOS2POSTGIS(geos_result, is3d);
3005 GEOSGeom_destroy(geos_result);
3006 if (!result)
3007 {
3008 elog(ERROR, "%s returned an error", __func__);
3009 PG_RETURN_NULL(); /*never get here */
3010 }
3011
3012 PG_RETURN_POINTER(result);
3013 }
3014
3015
3016 PG_FUNCTION_INFO_V1(clusterintersecting_garray);
clusterintersecting_garray(PG_FUNCTION_ARGS)3017 Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
3018 {
3019 Datum* result_array_data;
3020 ArrayType *array, *result;
3021 int is3d = 0;
3022 uint32 nelems, nclusters, i;
3023 GEOSGeometry **geos_inputs, **geos_results;
3024 int32_t srid = SRID_UNKNOWN;
3025
3026 /* Parameters used to construct a result array */
3027 int16 elmlen;
3028 bool elmbyval;
3029 char elmalign;
3030
3031 /* Null array, null geometry (should be empty?) */
3032 if (PG_ARGISNULL(0))
3033 PG_RETURN_NULL();
3034
3035 array = PG_GETARG_ARRAYTYPE_P(0);
3036 nelems = array_nelems_not_null(array);
3037
3038 POSTGIS_DEBUGF(3, "clusterintersecting_garray: number of non-null elements: %d", nelems);
3039
3040 if ( nelems == 0 ) PG_RETURN_NULL();
3041
3042 /* TODO short-circuit for one element? */
3043
3044 /* Ok, we really need geos now ;) */
3045 initGEOS(lwpgnotice, lwgeom_geos_error);
3046
3047 geos_inputs = ARRAY2GEOS(array, nelems, &is3d, &srid);
3048 if(!geos_inputs)
3049 {
3050 PG_RETURN_NULL();
3051 }
3052
3053 if (cluster_intersecting(geos_inputs, nelems, &geos_results, &nclusters) != LW_SUCCESS)
3054 {
3055 elog(ERROR, "clusterintersecting: Error performing clustering");
3056 PG_RETURN_NULL();
3057 }
3058 pfree(geos_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3059
3060 if (!geos_results) PG_RETURN_NULL();
3061
3062 result_array_data = palloc(nclusters * sizeof(Datum));
3063 for (i=0; i<nclusters; ++i)
3064 {
3065 result_array_data[i] = PointerGetDatum(GEOS2POSTGIS(geos_results[i], is3d));
3066 GEOSGeom_destroy(geos_results[i]);
3067 }
3068 lwfree(geos_results);
3069
3070 get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3071 result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3072
3073 if (!result)
3074 {
3075 elog(ERROR, "clusterintersecting: Error constructing return-array");
3076 PG_RETURN_NULL();
3077 }
3078
3079 PG_RETURN_POINTER(result);
3080 }
3081
3082 PG_FUNCTION_INFO_V1(cluster_within_distance_garray);
cluster_within_distance_garray(PG_FUNCTION_ARGS)3083 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS)
3084 {
3085 Datum* result_array_data;
3086 ArrayType *array, *result;
3087 int is3d = 0;
3088 uint32 nelems, nclusters, i;
3089 LWGEOM** lw_inputs;
3090 LWGEOM** lw_results;
3091 double tolerance;
3092 int32_t srid = SRID_UNKNOWN;
3093
3094 /* Parameters used to construct a result array */
3095 int16 elmlen;
3096 bool elmbyval;
3097 char elmalign;
3098
3099 /* Null array, null geometry (should be empty?) */
3100 if (PG_ARGISNULL(0))
3101 PG_RETURN_NULL();
3102
3103 array = PG_GETARG_ARRAYTYPE_P(0);
3104
3105 tolerance = PG_GETARG_FLOAT8(1);
3106 if (tolerance < 0)
3107 {
3108 lwpgerror("Tolerance must be a positive number.");
3109 PG_RETURN_NULL();
3110 }
3111
3112 nelems = array_nelems_not_null(array);
3113
3114 POSTGIS_DEBUGF(3, "cluster_within_distance_garray: number of non-null elements: %d", nelems);
3115
3116 if ( nelems == 0 ) PG_RETURN_NULL();
3117
3118 /* TODO short-circuit for one element? */
3119
3120 /* Ok, we really need geos now ;) */
3121 initGEOS(lwpgnotice, lwgeom_geos_error);
3122
3123 lw_inputs = ARRAY2LWGEOM(array, nelems, &is3d, &srid);
3124 if (!lw_inputs)
3125 {
3126 PG_RETURN_NULL();
3127 }
3128
3129 if (cluster_within_distance(lw_inputs, nelems, tolerance, &lw_results, &nclusters) != LW_SUCCESS)
3130 {
3131 elog(ERROR, "cluster_within: Error performing clustering");
3132 PG_RETURN_NULL();
3133 }
3134 pfree(lw_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3135
3136 if (!lw_results) PG_RETURN_NULL();
3137
3138 result_array_data = palloc(nclusters * sizeof(Datum));
3139 for (i=0; i<nclusters; ++i)
3140 {
3141 result_array_data[i] = PointerGetDatum(geometry_serialize(lw_results[i]));
3142 lwgeom_free(lw_results[i]);
3143 }
3144 lwfree(lw_results);
3145
3146 get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3147 result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3148
3149 if (!result)
3150 {
3151 elog(ERROR, "clusterwithin: Error constructing return-array");
3152 PG_RETURN_NULL();
3153 }
3154
3155 PG_RETURN_POINTER(result);
3156 }
3157
3158 PG_FUNCTION_INFO_V1(linemerge);
linemerge(PG_FUNCTION_ARGS)3159 Datum linemerge(PG_FUNCTION_ARGS)
3160 {
3161 GSERIALIZED *geom1;
3162 GSERIALIZED *result;
3163 LWGEOM *lwgeom1, *lwresult ;
3164
3165 geom1 = PG_GETARG_GSERIALIZED_P(0);
3166
3167
3168 lwgeom1 = lwgeom_from_gserialized(geom1) ;
3169
3170 lwresult = lwgeom_linemerge(lwgeom1);
3171 result = geometry_serialize(lwresult) ;
3172
3173 lwgeom_free(lwgeom1) ;
3174 lwgeom_free(lwresult) ;
3175
3176 PG_FREE_IF_COPY(geom1, 0);
3177
3178 PG_RETURN_POINTER(result);
3179 }
3180
3181 /*
3182 * Take a geometry and return an areal geometry
3183 * (Polygon or MultiPolygon).
3184 * Actually a wrapper around GEOSpolygonize,
3185 * transforming the resulting collection into
3186 * a valid polygon Geometry.
3187 */
3188 PG_FUNCTION_INFO_V1(ST_BuildArea);
ST_BuildArea(PG_FUNCTION_ARGS)3189 Datum ST_BuildArea(PG_FUNCTION_ARGS)
3190 {
3191 GSERIALIZED *result;
3192 GSERIALIZED *geom;
3193 LWGEOM *lwgeom_in, *lwgeom_out;
3194
3195 geom = PG_GETARG_GSERIALIZED_P(0);
3196 lwgeom_in = lwgeom_from_gserialized(geom);
3197
3198 lwgeom_out = lwgeom_buildarea(lwgeom_in);
3199 lwgeom_free(lwgeom_in) ;
3200
3201 if ( ! lwgeom_out )
3202 {
3203 PG_FREE_IF_COPY(geom, 0);
3204 PG_RETURN_NULL();
3205 }
3206
3207 result = geometry_serialize(lwgeom_out) ;
3208 lwgeom_free(lwgeom_out) ;
3209
3210 PG_FREE_IF_COPY(geom, 0);
3211 PG_RETURN_POINTER(result);
3212 }
3213
3214 /*
3215 * Take the vertices of a geometry and builds
3216 * Delaunay triangles around them.
3217 */
3218 PG_FUNCTION_INFO_V1(ST_DelaunayTriangles);
ST_DelaunayTriangles(PG_FUNCTION_ARGS)3219 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
3220 {
3221 GSERIALIZED *result;
3222 GSERIALIZED *geom;
3223 LWGEOM *lwgeom_in, *lwgeom_out;
3224 double tolerance = 0.0;
3225 int flags = 0;
3226
3227 geom = PG_GETARG_GSERIALIZED_P(0);
3228 tolerance = PG_GETARG_FLOAT8(1);
3229 flags = PG_GETARG_INT32(2);
3230
3231 lwgeom_in = lwgeom_from_gserialized(geom);
3232 lwgeom_out = lwgeom_delaunay_triangulation(lwgeom_in, tolerance, flags);
3233 lwgeom_free(lwgeom_in) ;
3234
3235 if ( ! lwgeom_out )
3236 {
3237 PG_FREE_IF_COPY(geom, 0);
3238 PG_RETURN_NULL();
3239 }
3240
3241 result = geometry_serialize(lwgeom_out) ;
3242 lwgeom_free(lwgeom_out) ;
3243
3244 PG_FREE_IF_COPY(geom, 0);
3245 PG_RETURN_POINTER(result);
3246 }
3247
3248 /*
3249 * ST_Snap
3250 *
3251 * Snap a geometry to another with a given tolerance
3252 */
3253 Datum ST_Snap(PG_FUNCTION_ARGS);
3254 PG_FUNCTION_INFO_V1(ST_Snap);
ST_Snap(PG_FUNCTION_ARGS)3255 Datum ST_Snap(PG_FUNCTION_ARGS)
3256 {
3257 GSERIALIZED *geom1, *geom2, *result;
3258 LWGEOM *lwgeom1, *lwgeom2, *lwresult;
3259 double tolerance;
3260
3261 geom1 = PG_GETARG_GSERIALIZED_P(0);
3262 geom2 = PG_GETARG_GSERIALIZED_P(1);
3263 tolerance = PG_GETARG_FLOAT8(2);
3264
3265 lwgeom1 = lwgeom_from_gserialized(geom1);
3266 lwgeom2 = lwgeom_from_gserialized(geom2);
3267
3268 lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3269 lwgeom_free(lwgeom1);
3270 lwgeom_free(lwgeom2);
3271 PG_FREE_IF_COPY(geom1, 0);
3272 PG_FREE_IF_COPY(geom2, 1);
3273
3274 result = geometry_serialize(lwresult);
3275 lwgeom_free(lwresult);
3276
3277 PG_RETURN_POINTER(result);
3278 }
3279
3280 /*
3281 * ST_Split
3282 *
3283 * Split polygon by line, line by line, line by point.
3284 * Returns at most components as a collection.
3285 * First element of the collection is always the part which
3286 * remains after the cut, while the second element is the
3287 * part which has been cut out. We arbitrarily take the part
3288 * on the *right* of cut lines as the part which has been cut out.
3289 * For a line cut by a point the part which remains is the one
3290 * from start of the line to the cut point.
3291 *
3292 *
3293 * Author: Sandro Santilli <strk@kbt.io>
3294 *
3295 * Work done for Faunalia (http://www.faunalia.it) with fundings
3296 * from Regione Toscana - Sistema Informativo per il Governo
3297 * del Territorio e dell'Ambiente (RT-SIGTA).
3298 *
3299 * Thanks to the PostGIS community for sharing poly/line ideas [1]
3300 *
3301 * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3302 *
3303 */
3304 Datum ST_Split(PG_FUNCTION_ARGS);
3305 PG_FUNCTION_INFO_V1(ST_Split);
ST_Split(PG_FUNCTION_ARGS)3306 Datum ST_Split(PG_FUNCTION_ARGS)
3307 {
3308 GSERIALIZED *in, *blade_in, *out;
3309 LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
3310
3311 in = PG_GETARG_GSERIALIZED_P(0);
3312 blade_in = PG_GETARG_GSERIALIZED_P(1);
3313 gserialized_error_if_srid_mismatch(in, blade_in, __func__);
3314
3315 lwgeom_in = lwgeom_from_gserialized(in);
3316 lwblade_in = lwgeom_from_gserialized(blade_in);
3317
3318 lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3319 lwgeom_free(lwgeom_in);
3320 lwgeom_free(lwblade_in);
3321
3322 if ( ! lwgeom_out )
3323 {
3324 PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3325 PG_FREE_IF_COPY(blade_in, 1);
3326 PG_RETURN_NULL();
3327 }
3328
3329 out = geometry_serialize(lwgeom_out);
3330 lwgeom_free(lwgeom_out);
3331 PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3332 PG_FREE_IF_COPY(blade_in, 1);
3333
3334 PG_RETURN_POINTER(out);
3335 }
3336
3337 /**********************************************************************
3338 *
3339 * ST_SharedPaths
3340 *
3341 * Return the set of paths shared between two linear geometries,
3342 * and their direction (same or opposite).
3343 *
3344 * Developed by Sandro Santilli (strk@kbt.io) for Faunalia
3345 * (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3346 * Informativo per la Gestione del Territorio e dell' Ambiente
3347 * [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3348 * trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3349 * 0494241492)"
3350 *
3351 **********************************************************************/
3352 Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3353 PG_FUNCTION_INFO_V1(ST_SharedPaths);
ST_SharedPaths(PG_FUNCTION_ARGS)3354 Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3355 {
3356 GSERIALIZED *geom1, *geom2, *out;
3357 LWGEOM *g1, *g2, *lwgeom_out;
3358
3359 geom1 = PG_GETARG_GSERIALIZED_P(0);
3360 geom2 = PG_GETARG_GSERIALIZED_P(1);
3361
3362 g1 = lwgeom_from_gserialized(geom1);
3363 g2 = lwgeom_from_gserialized(geom2);
3364
3365 lwgeom_out = lwgeom_sharedpaths(g1, g2);
3366 lwgeom_free(g1);
3367 lwgeom_free(g2);
3368
3369 if ( ! lwgeom_out )
3370 {
3371 PG_FREE_IF_COPY(geom1, 0);
3372 PG_FREE_IF_COPY(geom2, 1);
3373 PG_RETURN_NULL();
3374 }
3375
3376 out = geometry_serialize(lwgeom_out);
3377 lwgeom_free(lwgeom_out);
3378
3379 PG_FREE_IF_COPY(geom1, 0);
3380 PG_FREE_IF_COPY(geom2, 1);
3381 PG_RETURN_POINTER(out);
3382 }
3383
3384
3385 /**********************************************************************
3386 *
3387 * ST_Node
3388 *
3389 * Fully node a set of lines using the least possible nodes while
3390 * preserving all of the input ones.
3391 *
3392 **********************************************************************/
3393 Datum ST_Node(PG_FUNCTION_ARGS);
3394 PG_FUNCTION_INFO_V1(ST_Node);
ST_Node(PG_FUNCTION_ARGS)3395 Datum ST_Node(PG_FUNCTION_ARGS)
3396 {
3397 GSERIALIZED *geom1, *out;
3398 LWGEOM *g1, *lwgeom_out;
3399
3400 geom1 = PG_GETARG_GSERIALIZED_P(0);
3401
3402 g1 = lwgeom_from_gserialized(geom1);
3403
3404 lwgeom_out = lwgeom_node(g1);
3405 lwgeom_free(g1);
3406
3407 if ( ! lwgeom_out )
3408 {
3409 PG_FREE_IF_COPY(geom1, 0);
3410 PG_RETURN_NULL();
3411 }
3412
3413 out = geometry_serialize(lwgeom_out);
3414 lwgeom_free(lwgeom_out);
3415
3416 PG_FREE_IF_COPY(geom1, 0);
3417 PG_RETURN_POINTER(out);
3418 }
3419
3420 /******************************************
3421 *
3422 * ST_Voronoi
3423 *
3424 * Returns a Voronoi diagram constructed
3425 * from the points of the input geometry.
3426 *
3427 ******************************************/
3428 Datum ST_Voronoi(PG_FUNCTION_ARGS);
3429 PG_FUNCTION_INFO_V1(ST_Voronoi);
ST_Voronoi(PG_FUNCTION_ARGS)3430 Datum ST_Voronoi(PG_FUNCTION_ARGS)
3431 {
3432 GSERIALIZED* input;
3433 GSERIALIZED* clip;
3434 GSERIALIZED* result;
3435 LWGEOM* lwgeom_input;
3436 LWGEOM* lwgeom_result;
3437 double tolerance;
3438 GBOX clip_envelope;
3439 int custom_clip_envelope;
3440 int return_polygons;
3441
3442 /* Return NULL on NULL geometry */
3443 if (PG_ARGISNULL(0))
3444 PG_RETURN_NULL();
3445
3446 /* Read our tolerance value */
3447 if (PG_ARGISNULL(2))
3448 {
3449 lwpgerror("Tolerance must be a positive number.");
3450 PG_RETURN_NULL();
3451 }
3452
3453 tolerance = PG_GETARG_FLOAT8(2);
3454
3455 if (tolerance < 0)
3456 {
3457 lwpgerror("Tolerance must be a positive number.");
3458 PG_RETURN_NULL();
3459 }
3460
3461 /* Are we returning lines or polygons? */
3462 if (PG_ARGISNULL(3))
3463 {
3464 lwpgerror("return_polygons must be true or false.");
3465 PG_RETURN_NULL();
3466 }
3467 return_polygons = PG_GETARG_BOOL(3);
3468
3469 /* Read our clipping envelope, if applicable. */
3470 custom_clip_envelope = !PG_ARGISNULL(1);
3471 if (custom_clip_envelope) {
3472 clip = PG_GETARG_GSERIALIZED_P(1);
3473 if (!gserialized_get_gbox_p(clip, &clip_envelope))
3474 {
3475 lwpgerror("Could not determine envelope of clipping geometry.");
3476 PG_FREE_IF_COPY(clip, 1);
3477 PG_RETURN_NULL();
3478 }
3479 PG_FREE_IF_COPY(clip, 1);
3480 }
3481
3482 /* Read our input geometry */
3483 input = PG_GETARG_GSERIALIZED_P(0);
3484
3485 lwgeom_input = lwgeom_from_gserialized(input);
3486
3487 if(!lwgeom_input)
3488 {
3489 lwpgerror("Could not read input geometry.");
3490 PG_FREE_IF_COPY(input, 0);
3491 PG_RETURN_NULL();
3492 }
3493
3494 lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
3495 lwgeom_free(lwgeom_input);
3496
3497 if (!lwgeom_result)
3498 {
3499 lwpgerror("Error computing Voronoi diagram.");
3500 PG_FREE_IF_COPY(input, 0);
3501 PG_RETURN_NULL();
3502 }
3503
3504 result = geometry_serialize(lwgeom_result);
3505 lwgeom_free(lwgeom_result);
3506
3507 PG_FREE_IF_COPY(input, 0);
3508 PG_RETURN_POINTER(result);
3509 }
3510
3511 /******************************************
3512 *
3513 * ST_MinimumClearance
3514 *
3515 * Returns the minimum clearance of a geometry.
3516 *
3517 ******************************************/
3518 Datum ST_MinimumClearance(PG_FUNCTION_ARGS);
3519 PG_FUNCTION_INFO_V1(ST_MinimumClearance);
ST_MinimumClearance(PG_FUNCTION_ARGS)3520 Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
3521 {
3522 GSERIALIZED* input;
3523 GEOSGeometry* input_geos;
3524 int error;
3525 double result;
3526
3527 initGEOS(lwpgnotice, lwgeom_geos_error);
3528
3529 input = PG_GETARG_GSERIALIZED_P(0);
3530 input_geos = POSTGIS2GEOS(input);
3531 if (!input_geos)
3532 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3533
3534 error = GEOSMinimumClearance(input_geos, &result);
3535 GEOSGeom_destroy(input_geos);
3536 if (error) HANDLE_GEOS_ERROR("Error computing minimum clearance");
3537
3538 PG_FREE_IF_COPY(input, 0);
3539 PG_RETURN_FLOAT8(result);
3540 }
3541
3542 /******************************************
3543 *
3544 * ST_MinimumClearanceLine
3545 *
3546 * Returns the minimum clearance line of a geometry.
3547 *
3548 ******************************************/
3549 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS);
3550 PG_FUNCTION_INFO_V1(ST_MinimumClearanceLine);
ST_MinimumClearanceLine(PG_FUNCTION_ARGS)3551 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
3552 {
3553 GSERIALIZED* input;
3554 GSERIALIZED* result;
3555 GEOSGeometry* input_geos;
3556 GEOSGeometry* result_geos;
3557 int32_t srid;
3558
3559 initGEOS(lwpgnotice, lwgeom_geos_error);
3560
3561 input = PG_GETARG_GSERIALIZED_P(0);
3562 srid = gserialized_get_srid(input);
3563 input_geos = POSTGIS2GEOS(input);
3564 if (!input_geos)
3565 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3566
3567 result_geos = GEOSMinimumClearanceLine(input_geos);
3568 GEOSGeom_destroy(input_geos);
3569 if (!result_geos)
3570 HANDLE_GEOS_ERROR("Error computing minimum clearance");
3571
3572 GEOSSetSRID(result_geos, srid);
3573 result = GEOS2POSTGIS(result_geos, LW_FALSE);
3574 GEOSGeom_destroy(result_geos);
3575
3576 PG_FREE_IF_COPY(input, 0);
3577 PG_RETURN_POINTER(result);
3578 }
3579
3580 /******************************************
3581 *
3582 * ST_OrientedEnvelope
3583 *
3584 ******************************************/
3585 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS);
3586 PG_FUNCTION_INFO_V1(ST_OrientedEnvelope);
ST_OrientedEnvelope(PG_FUNCTION_ARGS)3587 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
3588 {
3589 GSERIALIZED* input;
3590 GSERIALIZED* result;
3591 GEOSGeometry* input_geos;
3592 GEOSGeometry* result_geos;
3593 int32_t srid;
3594
3595 initGEOS(lwpgnotice, lwgeom_geos_error);
3596
3597 input = PG_GETARG_GSERIALIZED_P(0);
3598 srid = gserialized_get_srid(input);
3599 input_geos = POSTGIS2GEOS(input);
3600 if (!input_geos)
3601 HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3602
3603 result_geos = GEOSMinimumRotatedRectangle(input_geos);
3604 GEOSGeom_destroy(input_geos);
3605 if (!result_geos)
3606 HANDLE_GEOS_ERROR("Error computing oriented envelope");
3607
3608 GEOSSetSRID(result_geos, srid);
3609 result = GEOS2POSTGIS(result_geos, LW_FALSE);
3610 GEOSGeom_destroy(result_geos);
3611
3612 PG_FREE_IF_COPY(input, 0);
3613 PG_RETURN_POINTER(result);
3614 }
3615