PostGIS  3.7.0dev-r@@SVN_REVISION@@
lwgeom_geos_predicates.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright 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_geos.h"
41 #include "liblwgeom.h"
42 #include "liblwgeom_internal.h"
43 #include "lwgeom_itree.h"
44 #include "lwgeom_geos_prepared.h"
45 #include "lwgeom_accum.h"
46 
47 
48 /* Prototypes for SQL-bound functions */
49 Datum relate_full(PG_FUNCTION_ARGS);
50 Datum relate_pattern(PG_FUNCTION_ARGS);
51 Datum disjoint(PG_FUNCTION_ARGS);
52 Datum touches(PG_FUNCTION_ARGS);
53 Datum ST_Intersects(PG_FUNCTION_ARGS);
54 Datum crosses(PG_FUNCTION_ARGS);
55 Datum contains(PG_FUNCTION_ARGS);
56 Datum within(PG_FUNCTION_ARGS);
57 Datum containsproperly(PG_FUNCTION_ARGS);
58 Datum covers(PG_FUNCTION_ARGS);
59 Datum overlaps(PG_FUNCTION_ARGS);
60 Datum coveredby(PG_FUNCTION_ARGS);
61 Datum ST_Equals(PG_FUNCTION_ARGS);
62 
63 
64 /*
65  * Utility to quickly check for polygonal geometries
66  */
67 static inline char
68 is_poly(const GSERIALIZED* g)
69 {
70  int type = gserialized_get_type(g);
71  return type == POLYGONTYPE || type == MULTIPOLYGONTYPE;
72 }
73 
74 /*
75  * Utility to quickly check for point geometries
76  */
77 static inline char
78 is_point(const GSERIALIZED* g)
79 {
80  int type = gserialized_get_type(g);
81  return type == POINTTYPE || type == MULTIPOINTTYPE;
82 }
83 
84 
85 
86 
87 
88 
89 
90 PG_FUNCTION_INFO_V1(ST_Intersects);
91 Datum ST_Intersects(PG_FUNCTION_ARGS)
92 {
93  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
94  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
95  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
96  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
97  int result;
98  GBOX box1, box2;
99  PrepGeomCache *prep_cache;
100 
101  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
102 
103  /* A.Intersects(Empty) == FALSE */
104  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
105  PG_RETURN_BOOL(false);
106 
107  /*
108  * Short-circuit 1: if geom2 bounding box does not overlap
109  * geom1 bounding box we can return FALSE.
110  */
111  if ( gserialized_get_gbox_p(geom1, &box1) &&
112  gserialized_get_gbox_p(geom2, &box2) )
113  {
114  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
115  PG_RETURN_BOOL(false);
116  }
117 
118  /*
119  * Short-circuit 2: if the geoms are a point and a polygon,
120  * call the itree_pip_intersects function.
121  */
122  if ((is_point(geom1) && is_poly(geom2)) ||
123  (is_point(geom2) && is_poly(geom1)))
124  {
125  SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
126  SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
127  const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
128  LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
129  IntervalTree *itree = GetIntervalTree(fcinfo, shared_gpoly);
130  bool result = itree_pip_intersects(itree, lwpt);
131  lwgeom_free(lwpt);
132  PG_RETURN_BOOL(result);
133  }
134 
135  initGEOS(lwpgnotice, lwgeom_geos_error);
136 
137  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
138  if ( prep_cache && prep_cache->prepared_geom )
139  {
140  GEOSGeometry *g = prep_cache->gcache.argnum == 1
141  ? POSTGIS2GEOS(geom2)
142  : POSTGIS2GEOS(geom1);
143  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
144  result = GEOSPreparedIntersects(prep_cache->prepared_geom, g);
145  GEOSGeom_destroy(g);
146  }
147  else
148  {
149  GEOSGeometry *g1, *g2;
150  g1 = POSTGIS2GEOS(geom1);
151  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
152  g2 = POSTGIS2GEOS(geom2);
153  if (!g2)
154  {
155  GEOSGeom_destroy(g1);
156  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
157  }
158  result = GEOSIntersects(g1, g2);
159  GEOSGeom_destroy(g1);
160  GEOSGeom_destroy(g2);
161  }
162 
163  if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
164 
165  PG_RETURN_BOOL(result);
166 }
167 
168 
169 PG_FUNCTION_INFO_V1(ST_Equals);
170 Datum ST_Equals(PG_FUNCTION_ARGS)
171 {
172  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
173  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
174  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
175  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
176  char result;
177  GBOX box1, box2;
178  GEOSGeometry *g1, *g2;
179 
180  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
181 
182  /* Empty == Empty */
183  if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
184  PG_RETURN_BOOL(true);
185 
186  /*
187  * Short-circuit: If geom1 and geom2 do not have the same bounding box
188  * we can return FALSE.
189  */
190  if ( gserialized_get_gbox_p(geom1, &box1) &&
191  gserialized_get_gbox_p(geom2, &box2) )
192  {
193  if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
194  {
195  PG_RETURN_BOOL(false);
196  }
197  }
198 
199  /*
200  * Short-circuit: if geom1 and geom2 are binary-equivalent, we can return
201  * TRUE. This is much faster than doing the comparison using GEOS.
202  */
203  if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
204  PG_RETURN_BOOL(true);
205  }
206 
207  initGEOS(lwpgnotice, lwgeom_geos_error);
208 
209  g1 = POSTGIS2GEOS(geom1);
210  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
211  g2 = POSTGIS2GEOS(geom2);
212  if (!g2)
213  {
214  GEOSGeom_destroy(g1);
215  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
216  }
217  result = GEOSEquals(g1, g2);
218  GEOSGeom_destroy(g1);
219  GEOSGeom_destroy(g2);
220 
221  if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
222 
223  PG_RETURN_BOOL(result);
224 }
225 
226 
227 
228 PG_FUNCTION_INFO_V1(touches);
229 Datum touches(PG_FUNCTION_ARGS)
230 {
231  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
232  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
233  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
234  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
235  char result;
236  GBOX box1, box2;
237  PrepGeomCache *prep_cache;
238 
239  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
240 
241  /* A.Touches(Empty) == FALSE */
242  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
243  PG_RETURN_BOOL(false);
244 
245  /*
246  * Short-circuit 1: if geom2 bounding box does not overlap
247  * geom1 bounding box we can return FALSE.
248  */
249  if ( gserialized_get_gbox_p(geom1, &box1) &&
250  gserialized_get_gbox_p(geom2, &box2) )
251  {
252  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
253  {
254  PG_RETURN_BOOL(false);
255  }
256  }
257 
258  initGEOS(lwpgnotice, lwgeom_geos_error);
259 
260  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
261  if ( prep_cache && prep_cache->prepared_geom )
262  {
263  GEOSGeometry *g = prep_cache->gcache.argnum == 1
264  ? POSTGIS2GEOS(geom2)
265  : POSTGIS2GEOS(geom1);
266  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
267  result = GEOSPreparedTouches(prep_cache->prepared_geom, g);
268  GEOSGeom_destroy(g);
269  }
270  else
271  {
272  GEOSGeometry *g1, *g2;
273  g1 = POSTGIS2GEOS(geom1);
274  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
275  g2 = POSTGIS2GEOS(geom2);
276  if (!g2)
277  {
278  GEOSGeom_destroy(g1);
279  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
280  }
281  result = GEOSTouches(g1, g2);
282  GEOSGeom_destroy(g1);
283  GEOSGeom_destroy(g2);
284  }
285 
286  if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
287 
288  PG_RETURN_BOOL(result);
289 }
290 
291 
292 PG_FUNCTION_INFO_V1(disjoint);
293 Datum disjoint(PG_FUNCTION_ARGS)
294 {
295  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
296  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
297  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
298  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
299  char result;
300  GBOX box1, box2;
301  PrepGeomCache *prep_cache;
302 
303  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
304 
305  /* A.Disjoint(Empty) == TRUE */
306  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
307  PG_RETURN_BOOL(true);
308 
309  /*
310  * Short-circuit 1: if geom2 bounding box does not overlap
311  * geom1 bounding box we can return TRUE.
312  */
313  if ( gserialized_get_gbox_p(geom1, &box1) &&
314  gserialized_get_gbox_p(geom2, &box2) )
315  {
316  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
317  {
318  PG_RETURN_BOOL(true);
319  }
320  }
321 
322  initGEOS(lwpgnotice, lwgeom_geos_error);
323 
324  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
325  if ( prep_cache && prep_cache->prepared_geom )
326  {
327  GEOSGeometry *g = prep_cache->gcache.argnum == 1
328  ? POSTGIS2GEOS(geom2)
329  : POSTGIS2GEOS(geom1);
330  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
331  result = GEOSPreparedDisjoint(prep_cache->prepared_geom, g);
332  GEOSGeom_destroy(g);
333  }
334  else
335  {
336  GEOSGeometry *g1, *g2;
337  g1 = POSTGIS2GEOS(geom1);
338  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
339  g2 = POSTGIS2GEOS(geom2);
340  if (!g2)
341  {
342  GEOSGeom_destroy(g1);
343  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
344  }
345  result = GEOSDisjoint(g1, g2);
346  GEOSGeom_destroy(g1);
347  GEOSGeom_destroy(g2);
348  }
349 
350  if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
351 
352  PG_RETURN_BOOL(result);
353 }
354 
355 
356 
360 PG_FUNCTION_INFO_V1(overlaps);
361 Datum overlaps(PG_FUNCTION_ARGS)
362 {
363  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
364  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
365  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
366  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
367  char result;
368  GBOX box1, box2;
369  PrepGeomCache *prep_cache;
370 
371  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
372 
373  /* A.Overlaps(Empty) == FALSE */
374  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
375  PG_RETURN_BOOL(false);
376 
377  /*
378  * Short-circuit 1: if geom2 bounding box does not overlap
379  * geom1 bounding box we can return FALSE.
380  */
381  if ( gserialized_get_gbox_p(geom1, &box1) &&
382  gserialized_get_gbox_p(geom2, &box2) )
383  {
384  if ( ! gbox_overlaps_2d(&box1, &box2) )
385  {
386  PG_RETURN_BOOL(false);
387  }
388  }
389 
390  initGEOS(lwpgnotice, lwgeom_geos_error);
391 
392  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
393  if ( prep_cache && prep_cache->prepared_geom )
394  {
395  GEOSGeometry *g = prep_cache->gcache.argnum == 1
396  ? POSTGIS2GEOS(geom2)
397  : POSTGIS2GEOS(geom1);
398  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
399  result = GEOSPreparedOverlaps(prep_cache->prepared_geom, g);
400  GEOSGeom_destroy(g);
401  }
402  else
403  {
404  GEOSGeometry *g1, *g2;
405  g1 = POSTGIS2GEOS(geom1);
406  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
407  g2 = POSTGIS2GEOS(geom2);
408  if (!g2)
409  {
410  GEOSGeom_destroy(g1);
411  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
412  }
413  result = GEOSOverlaps(g1, g2);
414  GEOSGeom_destroy(g1);
415  GEOSGeom_destroy(g2);
416  }
417 
418  if (result == 2) HANDLE_GEOS_ERROR("GEOSOverlaps");
419 
420  PG_RETURN_BOOL(result);
421 }
422 
423 
427 PG_FUNCTION_INFO_V1(crosses);
428 Datum crosses(PG_FUNCTION_ARGS)
429 {
430  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
431  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
432  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
433  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
434  int result;
435  GBOX box1, box2;
436  PrepGeomCache *prep_cache;
437 
438  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
439 
440  /* A.Crosses(Empty) == FALSE */
441  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
442  PG_RETURN_BOOL(false);
443 
444  /*
445  * Short-circuit 1: if geom2 bounding box does not overlap
446  * geom1 bounding box we can return FALSE.
447  */
448  if ( gserialized_get_gbox_p(geom1, &box1) &&
449  gserialized_get_gbox_p(geom2, &box2) )
450  {
451  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
452  {
453  PG_RETURN_BOOL(false);
454  }
455  }
456 
457  initGEOS(lwpgnotice, lwgeom_geos_error);
458 
459  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
460  if ( prep_cache && prep_cache->prepared_geom )
461  {
462  GEOSGeometry *g = prep_cache->gcache.argnum == 1
463  ? POSTGIS2GEOS(geom2)
464  : POSTGIS2GEOS(geom1);
465  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
466  result = GEOSPreparedCrosses(prep_cache->prepared_geom, g);
467  GEOSGeom_destroy(g);
468  }
469  else
470  {
471  GEOSGeometry *g1, *g2;
472  g1 = POSTGIS2GEOS(geom1);
473  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
474  g2 = POSTGIS2GEOS(geom2);
475  if (!g2)
476  {
477  GEOSGeom_destroy(g1);
478  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
479  }
480  result = GEOSCrosses(g1, g2);
481  GEOSGeom_destroy(g1);
482  GEOSGeom_destroy(g2);
483  }
484 
485  if (result == 2) HANDLE_GEOS_ERROR("GEOSCrosses");
486 
487  PG_RETURN_BOOL(result);
488 }
489 
490 
491 
495 PG_FUNCTION_INFO_V1(contains);
496 Datum contains(PG_FUNCTION_ARGS)
497 {
498  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
499  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
500  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
501  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
502  int result;
503  GEOSGeometry *g1, *g2;
504  GBOX box1, box2;
505  PrepGeomCache *prep_cache;
506 
507  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
508 
509  /* A.Contains(Empty) == FALSE */
510  if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
511  PG_RETURN_BOOL(false);
512 
513  POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
514 
515  /*
516  ** Short-circuit 1: if geom2 bounding box is not completely inside
517  ** geom1 bounding box we can return FALSE.
518  */
519  if (gserialized_get_gbox_p(geom1, &box1) &&
520  gserialized_get_gbox_p(geom2, &box2))
521  {
522  if (!gbox_contains_2d(&box1, &box2))
523  PG_RETURN_BOOL(false);
524  }
525 
526  /*
527  ** Short-circuit 2: if geom2 is a point and geom1 is a polygon
528  ** call the point-in-polygon function.
529  */
530  if (is_poly(geom1) && is_point(geom2))
531  {
532  const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom2);
533  LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
534  IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom1);
535  bool result = itree_pip_contains(itree, lwpt);
536  lwgeom_free(lwpt);
537  PG_RETURN_BOOL(result);
538  }
539 
540  initGEOS(lwpgnotice, lwgeom_geos_error);
541 
542  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, NULL);
543  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
544  {
545  g1 = POSTGIS2GEOS(geom2);
546  if (!g1)
547  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
548 
549  result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
550  GEOSGeom_destroy(g1);
551  }
552  else
553  {
554  g1 = POSTGIS2GEOS(geom1);
555  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
556  g2 = POSTGIS2GEOS(geom2);
557  if (!g2)
558  {
559  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
560  GEOSGeom_destroy(g1);
561  }
562  result = GEOSContains( g1, g2);
563  GEOSGeom_destroy(g1);
564  GEOSGeom_destroy(g2);
565  }
566 
567  if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
568 
569  PG_RETURN_BOOL(result > 0);
570 }
571 
572 
576 PG_FUNCTION_INFO_V1(within);
577 Datum within(PG_FUNCTION_ARGS)
578 {
579  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
580  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
581  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
582  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
583  int result;
584  GEOSGeometry *g1, *g2;
585  GBOX box1, box2;
586  PrepGeomCache *prep_cache;
587 
588  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
589 
590  /* A.Contains(Empty) == FALSE */
591  if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
592  PG_RETURN_BOOL(false);
593 
594  POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
595 
596  /*
597  ** Short-circuit 1: if geom1 bounding box is not completely inside
598  ** geom2 bounding box we can return FALSE.
599  */
600  if (gserialized_get_gbox_p(geom1, &box1) &&
601  gserialized_get_gbox_p(geom2, &box2))
602  {
603  if (!gbox_within_2d(&box1, &box2))
604  PG_RETURN_BOOL(false);
605  }
606 
607  /*
608  ** Short-circuit 2: if geom2 is a polygon and geom1 is a point
609  ** call the point-in-polygon function.
610  */
611  if (is_poly(geom2) && is_point(geom1))
612  {
613  const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom1);
614  LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
615  IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom2);
616  bool result = itree_pip_contains(itree, lwpt);
617  lwgeom_free(lwpt);
618  PG_RETURN_BOOL(result);
619  }
620 
621  initGEOS(lwpgnotice, lwgeom_geos_error);
622 
623  prep_cache = GetPrepGeomCache(fcinfo, NULL, shared_geom2);
624  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 2 )
625  {
626  g1 = POSTGIS2GEOS(geom1);
627  if (!g1)
628  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
629 
630  result = GEOSPreparedContains(prep_cache->prepared_geom, g1);
631  GEOSGeom_destroy(g1);
632  }
633  else
634  {
635  g1 = POSTGIS2GEOS(geom1);
636  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
637  g2 = POSTGIS2GEOS(geom2);
638  if (!g2)
639  {
640  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
641  GEOSGeom_destroy(g1);
642  }
643  result = GEOSWithin(g1, g2);
644  GEOSGeom_destroy(g1);
645  GEOSGeom_destroy(g2);
646  }
647 
648  if ( result == 2 ) HANDLE_GEOS_ERROR("GEOSWithin");
649 
650  PG_RETURN_BOOL(result > 0);
651 }
652 
656 PG_FUNCTION_INFO_V1(containsproperly);
657 Datum containsproperly(PG_FUNCTION_ARGS)
658 {
659  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
660  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
661  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
662  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
663  char result;
664  GBOX box1, box2;
665  PrepGeomCache *prep_cache;
666 
667  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
668 
669  /* A.ContainsProperly(Empty) == FALSE */
670  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
671  PG_RETURN_BOOL(false);
672 
673  /*
674  * Short-circuit: if geom2 bounding box is not completely inside
675  * geom1 bounding box we can return FALSE.
676  */
677  if ( gserialized_get_gbox_p(geom1, &box1) &&
678  gserialized_get_gbox_p(geom2, &box2) )
679  {
680  if ( ! gbox_contains_2d(&box1, &box2) )
681  PG_RETURN_BOOL(false);
682  }
683 
684  initGEOS(lwpgnotice, lwgeom_geos_error);
685 
686  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
687  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
688  {
689  GEOSGeometry *g = POSTGIS2GEOS(geom2);
690  if (!g) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
691  result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
692  GEOSGeom_destroy(g);
693  }
694  else
695  {
696  GEOSGeometry *g1, *g2;
697  g1 = POSTGIS2GEOS(geom1);
698  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
699  g2 = POSTGIS2GEOS(geom2);
700  if (!g2)
701  {
702  GEOSGeom_destroy(g1);
703  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
704  }
705  result = GEOSRelatePattern(g1, g2, "T**FF*FF*");
706  GEOSGeom_destroy(g1);
707  GEOSGeom_destroy(g2);
708  }
709 
710  if (result == 2) HANDLE_GEOS_ERROR("GEOSContainProperly");
711 
712  PG_RETURN_BOOL(result);
713 }
714 
715 
721 PG_FUNCTION_INFO_V1(covers);
722 Datum covers(PG_FUNCTION_ARGS)
723 {
724  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
725  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
726  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
727  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
728  int result;
729  GBOX box1, box2;
730  PrepGeomCache *prep_cache;
731 
732  POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
733 
734  /* A.Covers(Empty) == FALSE */
735  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
736  PG_RETURN_BOOL(false);
737 
738  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
739 
740  /*
741  * Short-circuit 1: if geom2 bounding box is not completely inside
742  * geom1 bounding box we can return FALSE.
743  */
744  if ( gserialized_get_gbox_p(geom1, &box1) &&
745  gserialized_get_gbox_p(geom2, &box2) )
746  {
747  if ( ! gbox_contains_2d(&box1, &box2) )
748  {
749  PG_RETURN_BOOL(false);
750  }
751  }
752  /*
753  * Short-circuit 2: if geom2 is a point and geom1 is a polygon
754  * call the point-in-polygon function.
755  */
756  if (is_poly(geom1) && is_point(geom2))
757  {
758  const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom2);
759  LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
760  IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom1);
761  bool result = itree_pip_covers(itree, lwpt);
762  lwgeom_free(lwpt);
763  PG_RETURN_BOOL(result);
764  }
765 
766  initGEOS(lwpgnotice, lwgeom_geos_error);
767 
768  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
769  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
770  {
771  GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
772  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
773  result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
774  GEOSGeom_destroy(g1);
775  }
776  else
777  {
778  GEOSGeometry *g1, *g2;
779 
780  g1 = POSTGIS2GEOS(geom1);
781  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
782  g2 = POSTGIS2GEOS(geom2);
783  if (!g2)
784  {
785  GEOSGeom_destroy(g1);
786  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
787  }
788  result = GEOSRelatePattern( g1, g2, "******FF*" );
789  GEOSGeom_destroy(g1);
790  GEOSGeom_destroy(g2);
791  }
792 
793  if (result == 2) HANDLE_GEOS_ERROR("GEOSCovers");
794 
795  PG_RETURN_BOOL(result);
796 }
797 
798 
805 PG_FUNCTION_INFO_V1(coveredby);
806 Datum coveredby(PG_FUNCTION_ARGS)
807 {
808  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
809  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
810  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
811  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
812  int result;
813  GBOX box1, box2;
814  PrepGeomCache *prep_cache;
815 
816  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
817 
818  POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
819 
820  /* A.CoveredBy(Empty) == FALSE */
821  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
822  PG_RETURN_BOOL(false);
823 
824  /*
825  * Short-circuit 1: if geom1 bounding box is not completely inside
826  * geom2 bounding box we can return FALSE.
827  */
828  if ( gserialized_get_gbox_p(geom1, &box1) &&
829  gserialized_get_gbox_p(geom2, &box2) )
830  {
831  if ( ! gbox_within_2d(&box1, &box2) )
832  {
833  PG_RETURN_BOOL(false);
834  }
835  }
836  /*
837  * Short-circuit 2: if geom1 is a point and geom2 is a polygon
838  * call the point-in-polygon function.
839  */
840  if (is_point(geom1) && is_poly(geom2))
841  {
842  const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom1);
843  IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom2);
844  LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
845  bool result = itree_pip_covers(itree, lwpt);
846  lwgeom_free(lwpt);
847  PG_RETURN_BOOL(result);
848  }
849 
850  initGEOS(lwpgnotice, lwgeom_geos_error);
851 
852  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
853  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 2)
854  {
855  GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
856  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
857  result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
858  GEOSGeom_destroy(g1);
859  }
860  else
861  {
862  GEOSGeometry *g1, *g2;
863  g1 = POSTGIS2GEOS(geom1);
864  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
865 
866  g2 = POSTGIS2GEOS(geom2);
867  if (!g2)
868  {
869  GEOSGeom_destroy(g1);
870  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
871  }
872 
873  result = GEOSCoveredBy(g1, g2);
874  GEOSGeom_destroy(g1);
875  GEOSGeom_destroy(g2);
876  }
877 
878  if (result == 2) HANDLE_GEOS_ERROR("GEOSCoveredBy");
879 
880  PG_RETURN_BOOL(result);
881 }
882 
883 #if POSTGIS_GEOS_VERSION >= 31300
884 /*
885  * Flip the DE9IM matrix around the diagonal to
886  * account for flipping the order of the operands
887  * to the GEOSPreparedRelatePattern function.
888  */
889 static void
890 imInvert(char *im)
891 {
892  char t;
893  t = im[1]; im[1] = im[3]; im[3] = t;
894  t = im[2]; im[2] = im[6]; im[6] = t;
895  t = im[5]; im[5] = im[7]; im[7] = t;
896 }
897 #endif
898 
899 PG_FUNCTION_INFO_V1(relate_pattern);
900 Datum relate_pattern(PG_FUNCTION_ARGS)
901 {
902  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
903  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
904  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
905  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
906 
907  /* Ensure DE9IM pattern is no more than 9 chars */
908  text *imPtr = DatumGetTextP(DirectFunctionCall2(text_left,
909  PG_GETARG_DATUM(2), Int32GetDatum(9)));
910  char *im = text_to_cstring(imPtr);
911  char result;
912  uint32_t i;
913 #if POSTGIS_GEOS_VERSION >= 31300
914  PrepGeomCache *prep_cache;
915 #endif
916 
917  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
918 
919  /*
920  ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
921  */
922  for (i = 0; i < strlen(im); i++)
923  im[i] = toupper(im[i]);
924 
925  initGEOS(lwpgnotice, lwgeom_geos_error);
926 
927 #if POSTGIS_GEOS_VERSION >= 31300
928  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
929  if ( prep_cache && prep_cache->prepared_geom )
930  {
931  GEOSGeometry *g = NULL;
932  if (prep_cache->gcache.argnum == 1)
933  {
934  g = POSTGIS2GEOS(geom2);
935  }
936  else
937  {
938  g = POSTGIS2GEOS(geom1);
939  imInvert(im);
940  }
941 
942  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
943  result = GEOSPreparedRelatePattern(prep_cache->prepared_geom, g, im);
944  GEOSGeom_destroy(g);
945  }
946  else
947 #endif
948  {
949  GEOSGeometry *g1, *g2;
950  g1 = POSTGIS2GEOS(geom1);
951  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
952  g2 = POSTGIS2GEOS(geom2);
953  if (!g2)
954  {
955  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
956  GEOSGeom_destroy(g1);
957  }
958  result = GEOSRelatePattern(g1, g2, im);
959  GEOSGeom_destroy(g1);
960  GEOSGeom_destroy(g2);
961  }
962 
963  pfree(im);
964  if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
965 
966  PG_RETURN_BOOL(result);
967 }
968 
969 
970 PG_FUNCTION_INFO_V1(relate_full);
971 Datum relate_full(PG_FUNCTION_ARGS)
972 {
973  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
974  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
975  GEOSGeometry *g1, *g2;
976  char *relate_str;
977  text *result;
978  int bnr = GEOSRELATE_BNR_OGC;
979 
980  /* TODO handle empty */
981  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
982 
983  if ( PG_NARGS() > 2 )
984  bnr = PG_GETARG_INT32(2);
985 
986  initGEOS(lwpgnotice, lwgeom_geos_error);
987 
988  g1 = POSTGIS2GEOS(geom1);
989  if (!g1)
990  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
991  g2 = POSTGIS2GEOS(geom2);
992  if (!g2)
993  {
994  GEOSGeom_destroy(g1);
995  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
996  }
997 
998  POSTGIS_DEBUG(3, "constructed geometries ");
999 
1000  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
1001  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
1002 
1003  relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
1004  if (!relate_str) {
1005  GEOSGeom_destroy(g1);
1006  GEOSGeom_destroy(g2);
1007  HANDLE_GEOS_ERROR("GEOSRelate");
1008  }
1009  result = cstring_to_text(relate_str);
1010  GEOSFree(relate_str);
1011  GEOSGeom_destroy(g1);
1012  GEOSGeom_destroy(g2);
1013 
1014  PG_FREE_IF_COPY(geom1, 0);
1015  PG_FREE_IF_COPY(geom2, 1);
1016 
1017  PG_RETURN_TEXT_P(result);
1018 }
1019 
1020 
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
gbox_overlaps_2d
int gbox_overlaps_2d(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the GBOX overlaps on the 2d plane, LW_FALSE otherwise.
Definition: gbox.c:323
gbox_within_2d
int gbox_within_2d(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the first GBOX is within the second on the 2d plane, LW_FALSE otherwise.
Definition: gbox.c:351
gbox_same_2d_float
int gbox_same_2d_float(const GBOX *g1, const GBOX *g2)
Check if two given GBOX are the same in x and y, or would round to the same GBOX in x and if serializ...
Definition: gbox.c:187
gbox_contains_2d
int gbox_contains_2d(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the first GBOX contains the second on the 2d plane, LW_FALSE otherwise.
Definition: gbox.c:339
gserialized_error_if_srid_mismatch
void gserialized_error_if_srid_mismatch(const GSERIALIZED *g1, const GSERIALIZED *g2, const char *funcname)
Definition: gserialized.c:432
gserialized_get_gbox_p
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *gbox)
Read the box from the GSERIALIZED or calculate it if necessary.
Definition: gserialized.c:94
lwgeom_from_gserialized
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:268
gserialized_is_empty
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: gserialized.c:181
gserialized_get_type
uint32_t gserialized_get_type(const GSERIALIZED *g)
Extract the geometry type from the serialized form (it hides in the anonymous data area,...
Definition: gserialized.c:118
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition: liblwgeom/lwgeom_geos.c:57
lwgeom_geos.h
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:94
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1218
MULTIPOINTTYPE
#define MULTIPOINTTYPE
Definition: liblwgeom.h:105
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:102
MULTIPOLYGONTYPE
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:107
POLYGONTYPE
#define POLYGONTYPE
Definition: liblwgeom.h:104
liblwgeom.h
This library is the generic geometry handling section of PostGIS.
liblwgeom_internal.h
lwgeom_accum.h
contains
Datum contains(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:496
coveredby
Datum coveredby(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:806
is_point
static char is_point(const GSERIALIZED *g)
Definition: lwgeom_geos_predicates.c:78
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:722
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(ST_Intersects)
touches
Datum touches(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:229
disjoint
Datum disjoint(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:293
ST_Intersects
Datum ST_Intersects(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:91
crosses
Datum crosses(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:428
relate_full
Datum relate_full(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:971
relate_pattern
Datum relate_pattern(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:900
is_poly
static char is_poly(const GSERIALIZED *g)
Definition: lwgeom_geos_predicates.c:68
overlaps
Datum overlaps(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:361
ST_Equals
Datum ST_Equals(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:170
within
Datum within(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:577
containsproperly
Datum containsproperly(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:657
GetPrepGeomCache
PrepGeomCache * GetPrepGeomCache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1, SHARED_GSERIALIZED *g2)
Given a couple potential geometries and a function call context, return a prepared structure for one ...
Definition: lwgeom_geos_prepared.c:369
lwgeom_geos_prepared.h
itree_pip_intersects
bool itree_pip_intersects(const IntervalTree *itree, const LWGEOM *lwpoints)
Definition: lwgeom_itree.c:217
GetIntervalTree
IntervalTree * GetIntervalTree(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1)
Checks for a cache hit against the provided geometry and returns a pre-built index structure (RTREE_P...
Definition: lwgeom_itree.c:114
itree_pip_covers
bool itree_pip_covers(const IntervalTree *itree, const LWGEOM *lwpoints)
Definition: lwgeom_itree.c:185
itree_pip_contains
bool itree_pip_contains(const IntervalTree *itree, const LWGEOM *lwpoints)
Definition: lwgeom_itree.c:140
lwgeom_itree.h
ovdump.type
type
Definition: ovdump.py:42
POSTGIS2GEOS
GEOSGeometry * POSTGIS2GEOS(const GSERIALIZED *pglwgeom)
Definition: postgis/lwgeom_geos.c:1822
HANDLE_GEOS_ERROR
#define HANDLE_GEOS_ERROR(label)
Definition: postgis/lwgeom_geos.h:56
PrepGeomCache::prepared_geom
const GEOSPreparedGeometry * prepared_geom
Definition: lwgeom_geos_prepared.h:54