PostGIS  3.7.0dev-r@@SVN_REVISION@@
gserialized_gist_2d.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 Paul Ramsey <pramsey@cleverelephant.ca>
22  * Copyright 2017-2025 Darafei Praliaskouski <me@komzpa.net>
23  *
24  **********************************************************************/
25 
26 /*
27 ** R-Tree Bibliography
28 **
29 ** [1] A. Guttman. R-tree: a dynamic index structure for spatial searching.
30 ** Proceedings of the ACM SIGMOD Conference, pp 47-57, June 1984.
31 ** [2] C.-H. Ang and T. C. Tan. New linear node splitting algorithm for
32 ** R-Trees. Advances in Spatial Databases - 5th International Symposium,
33 ** 1997
34 ** [3] N. Beckmann, H.-P. Kriegel, R. Schneider, B. Seeger. The R*tree: an
35 ** efficient and robust access method for points and rectangles.
36 ** Proceedings of the ACM SIGMOD Conference. June 1990.
37 ** [4] A. Korotkov, "A new double sorting-based node splitting algorithm for R-tree",
38 ** http://syrcose.ispras.ru/2011/files/SYRCoSE2011_Proceedings.pdf#page=36
39 */
40 
41 #include "postgres.h"
42 #include "access/gist.h" /* For GiST */
43 #include "access/itup.h"
44 #include "access/skey.h"
45 #include "utils/sortsupport.h" /* For index building sort support */
46 
47 #include "../postgis_config.h"
48 
49 #include "liblwgeom.h" /* For standard geometry types. */
50 #include "liblwgeom_internal.h"
51 #include "lwgeom_pg.h" /* For debugging macros. */
52 #include "gserialized_gist.h" /* For utility functions. */
53 
54 #include <float.h> /* For FLT_MAX */
55 #include <math.h>
56 
57 /* When is a node split not so good? Prefer to not split 2 pages into more than 3 pages in index. */
58 #define LIMIT_RATIO 0.3333333333333333
59 
60 /* How many index tuples does one page fit? Recursive splits will target this. */
61 #define INDEX_TUPLES_PER_PAGE 291
62 
63 /*
64 ** For debugging
65 */
66 #if POSTGIS_DEBUG_LEVEL > 0
67 static int g2d_counter_leaf = 0;
68 static int g2d_counter_internal = 0;
69 #endif
70 
71 /*
72 ** GiST 2D key stubs
73 */
74 Datum box2df_out(PG_FUNCTION_ARGS);
75 Datum box2df_in(PG_FUNCTION_ARGS);
76 
77 /*
78 ** GiST 2D index function prototypes
79 */
80 Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS);
81 Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS);
82 Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS);
83 Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS);
84 Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS);
85 Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS);
86 Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS);
87 Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS);
88 Datum gserialized_gist_sortsupport_2d(PG_FUNCTION_ARGS);
89 
90 /*
91 ** GiST 2D operator prototypes
92 */
93 Datum gserialized_same_2d(PG_FUNCTION_ARGS);
94 Datum gserialized_within_2d(PG_FUNCTION_ARGS);
95 Datum gserialized_contains_2d(PG_FUNCTION_ARGS);
96 Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS);
97 Datum gserialized_left_2d(PG_FUNCTION_ARGS);
98 Datum gserialized_right_2d(PG_FUNCTION_ARGS);
99 Datum gserialized_above_2d(PG_FUNCTION_ARGS);
100 Datum gserialized_below_2d(PG_FUNCTION_ARGS);
101 Datum gserialized_overleft_2d(PG_FUNCTION_ARGS);
102 Datum gserialized_overright_2d(PG_FUNCTION_ARGS);
103 Datum gserialized_overabove_2d(PG_FUNCTION_ARGS);
104 Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS);
105 Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS);
106 Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS);
107 Datum gserialized_contains_box2df_geom_2d(PG_FUNCTION_ARGS);
108 Datum gserialized_contains_box2df_box2df_2d(PG_FUNCTION_ARGS);
109 Datum gserialized_within_box2df_geom_2d(PG_FUNCTION_ARGS);
110 Datum gserialized_within_box2df_box2df_2d(PG_FUNCTION_ARGS);
111 Datum gserialized_overlaps_box2df_geom_2d(PG_FUNCTION_ARGS);
112 Datum gserialized_overlaps_box2df_box2df_2d(PG_FUNCTION_ARGS);
113 
114 /*
115 ** true/false test function type
116 */
117 typedef bool (*box2df_predicate)(const BOX2DF *a, const BOX2DF *b);
118 
119 
120 static char* box2df_to_string(const BOX2DF *a)
121 {
122  static const int precision = 12;
123  char tmp[13 + 4 * OUT_MAX_BYTES_DOUBLE] = {'B', 'O', 'X', '2', 'D', 'F', '(', 0};
124  int len = 7;
125 
126  if (a == NULL)
127  return pstrdup("<NULLPTR>");
128 
129  len += lwprint_double(a->xmin, precision, &tmp[len]);
130  tmp[len++] = ' ';
131  len += lwprint_double(a->ymin, precision, &tmp[len]);
132  tmp[len++] = ',';
133  tmp[len++] = ' ';
134  len += lwprint_double(a->xmax, precision, &tmp[len]);
135  tmp[len++] = ' ';
136  len += lwprint_double(a->ymax, precision, &tmp[len]);
137  tmp[len++] = ')';
138 
139  return pstrdup(tmp);
140 }
141 
142 /* Allocate a new copy of BOX2DF */
143 BOX2DF* box2df_copy(BOX2DF *b)
144 {
145  BOX2DF *c = (BOX2DF*)palloc(sizeof(BOX2DF));
146  memcpy((void*)c, (void*)b, sizeof(BOX2DF));
147  POSTGIS_DEBUGF(5, "copied box2df (%p) to box2df (%p)", b, c);
148  return c;
149 }
150 
151 inline bool box2df_is_empty(const BOX2DF *a)
152 {
153  if (isnan(a->xmin))
154  return true;
155  else
156  return false;
157 }
158 
159 inline void box2df_set_empty(BOX2DF *a)
160 {
161  a->xmin = a->xmax = a->ymin = a->ymax = NAN;
162  return;
163 }
164 
165 inline void box2df_set_finite(BOX2DF *a)
166 {
167  if ( ! isfinite(a->xmax) )
168  a->xmax = FLT_MAX;
169  if ( ! isfinite(a->ymax) )
170  a->ymax = FLT_MAX;
171  if ( ! isfinite(a->ymin) )
172  a->ymin = -1*FLT_MAX;
173  if ( ! isfinite(a->xmin) )
174  a->xmin = -1*FLT_MAX;
175  return;
176 }
177 
178 /* Enlarge b_union to contain b_new. */
179 void box2df_merge(BOX2DF *b_union, BOX2DF *b_new)
180 {
181 
182  POSTGIS_DEBUGF(5, "merging %s with %s", box2df_to_string(b_union), box2df_to_string(b_new));
183  /* Adjust minimums */
184  if (b_union->xmin > b_new->xmin || isnan(b_union->xmin))
185  b_union->xmin = b_new->xmin;
186  if (b_union->ymin > b_new->ymin || isnan(b_union->ymin))
187  b_union->ymin = b_new->ymin;
188  /* Adjust maximums */
189  if (b_union->xmax < b_new->xmax || isnan(b_union->xmax))
190  b_union->xmax = b_new->xmax;
191  if (b_union->ymax < b_new->ymax || isnan(b_union->ymax))
192  b_union->ymax = b_new->ymax;
193 
194  POSTGIS_DEBUGF(5, "merge complete %s", box2df_to_string(b_union));
195  return;
196 }
197 
198 
199 /* Convert a double-based GBOX into a float-based BOX2DF,
200  ensuring the float box is larger than the double box */
201 static inline int box2df_from_gbox_p(GBOX *box, BOX2DF *a)
202 {
203  a->xmin = next_float_down(box->xmin);
204  a->xmax = next_float_up(box->xmax);
205  a->ymin = next_float_down(box->ymin);
206  a->ymax = next_float_up(box->ymax);
207  return LW_SUCCESS;
208 }
209 
210 int box2df_to_gbox_p(BOX2DF *a, GBOX *box)
211 {
212  memset(box, 0, sizeof(GBOX));
213  box->xmin = a->xmin;
214  box->xmax = a->xmax;
215  box->ymin = a->ymin;
216  box->ymax = a->ymax;
217  return LW_SUCCESS;
218 }
219 
220 /***********************************************************************
221 ** BOX2DF tests for 2D index operators.
222 */
223 
224 /* Ensure all minimums are below maximums. */
225 inline void box2df_validate(BOX2DF *b)
226 {
227  float tmp;
228 
229  if ( box2df_is_empty(b) )
230  return;
231 
232  if ( b->xmax < b->xmin )
233  {
234  tmp = b->xmin;
235  b->xmin = b->xmax;
236  b->xmax = tmp;
237  }
238  if ( b->ymax < b->ymin )
239  {
240  tmp = b->ymin;
241  b->ymin = b->ymax;
242  b->ymax = tmp;
243  }
244  return;
245 }
246 
247 bool box2df_overlaps(const BOX2DF *a, const BOX2DF *b)
248 {
249  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
250  return false;
251 
252  if ( (a->xmin > b->xmax) || (b->xmin > a->xmax) ||
253  (a->ymin > b->ymax) || (b->ymin > a->ymax) )
254  {
255  return false;
256  }
257 
258  return true;
259 }
260 
261 bool box2df_contains(const BOX2DF *a, const BOX2DF *b)
262 {
263  if ( !a || !b )
264  return false;
265 
266  /* All things can contain EMPTY (except EMPTY) */
267  if ( box2df_is_empty(b) && ! box2df_is_empty(a) )
268  return true;
269 
270  if ( (a->xmin > b->xmin) || (a->xmax < b->xmax) ||
271  (a->ymin > b->ymin) || (a->ymax < b->ymax) )
272  {
273  return false;
274  }
275 
276  return true;
277 }
278 
279 static bool box2df_within(const BOX2DF *a, const BOX2DF *b)
280 {
281  if ( !a || !b )
282  return false;
283 
284  /* EMPTY is within all other things (except EMPTY) */
285  if ( box2df_is_empty(a) && ! box2df_is_empty(b) )
286  return true;
287 
288  if ( (a->xmin < b->xmin) || (a->xmax > b->xmax) ||
289  (a->ymin < b->ymin) || (a->ymax > b->ymax) )
290  {
291  return false;
292  }
293 
294  return true;
295 }
296 
297 bool box2df_equals(const BOX2DF *a, const BOX2DF *b)
298 {
299  if ( !a && !b )
300  return true;
301  else if ( !a || !b )
302  return false;
303  else if ( box2df_is_empty(a) && box2df_is_empty(b) )
304  return true;
305  else if ( box2df_is_empty(a) || box2df_is_empty(b) )
306  return false;
307  else if ((a->xmin == b->xmin) && (a->xmax == b->xmax) && (a->ymin == b->ymin) && (a->ymax == b->ymax))
308  return true;
309  else
310  return false;
311 }
312 
313 bool box2df_overleft(const BOX2DF *a, const BOX2DF *b)
314 {
315  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
316  return false;
317 
318  /* a.xmax <= b.xmax */
319  return a->xmax <= b->xmax;
320 }
321 
322 bool box2df_left(const BOX2DF *a, const BOX2DF *b)
323 {
324  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
325  return false;
326 
327  /* a.xmax < b.xmin */
328  return a->xmax < b->xmin;
329 }
330 
331 bool box2df_right(const BOX2DF *a, const BOX2DF *b)
332 {
333  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
334  return false;
335 
336  /* a.xmin > b.xmax */
337  return a->xmin > b->xmax;
338 }
339 
340 bool box2df_overright(const BOX2DF *a, const BOX2DF *b)
341 {
342  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
343  return false;
344 
345  /* a.xmin >= b.xmin */
346  return a->xmin >= b->xmin;
347 }
348 
349 bool box2df_overbelow(const BOX2DF *a, const BOX2DF *b)
350 {
351  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
352  return false;
353 
354  /* a.ymax <= b.ymax */
355  return a->ymax <= b->ymax;
356 }
357 
358 bool box2df_below(const BOX2DF *a, const BOX2DF *b)
359 {
360  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
361  return false;
362 
363  /* a.ymax < b.ymin */
364  return a->ymax < b->ymin;
365 }
366 
367 bool box2df_above(const BOX2DF *a, const BOX2DF *b)
368 {
369  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
370  return false;
371 
372  /* a.ymin > b.ymax */
373  return a->ymin > b->ymax;
374 }
375 
376 bool box2df_overabove(const BOX2DF *a, const BOX2DF *b)
377 {
378  if ( !a || !b || box2df_is_empty(a) || box2df_is_empty(b) )
379  return false;
380 
381  /* a.ymin >= b.ymin */
382  return a->ymin >= b->ymin;
383 }
384 
388 static double box2df_distance_leaf_centroid(const BOX2DF *a, const BOX2DF *b)
389 {
390  /* The centroid->centroid distance between the boxes */
391  double a_x = (a->xmax + a->xmin) / 2.0;
392  double a_y = (a->ymax + a->ymin) / 2.0;
393  double b_x = (b->xmax + b->xmin) / 2.0;
394  double b_y = (b->ymax + b->ymin) / 2.0;
395 
396  /* This "distance" is only used for comparisons, */
397  return sqrt((a_x - b_x) * (a_x - b_x) + (a_y - b_y) * (a_y - b_y));
398 }
399 
400 /* Quick distance function */
401 static inline double pt_distance(double ax, double ay, double bx, double by)
402 {
403  return sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
404 }
405 
409 static double box2df_distance(const BOX2DF *a, const BOX2DF *b)
410 {
411  /* Check for overlap */
412  if ( box2df_overlaps(a, b) )
413  return 0.0;
414 
415  if ( box2df_left(a, b) )
416  {
417  if ( box2df_above(a, b) )
418  return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
419  if ( box2df_below(a, b) )
420  return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
421  else
422  return (double)b->xmin - (double)a->xmax;
423  }
424  if ( box2df_right(a, b) )
425  {
426  if ( box2df_above(a, b) )
427  return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
428  if ( box2df_below(a, b) )
429  return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
430  else
431  return (double)a->xmin - (double)b->xmax;
432  }
433  if ( box2df_above(a, b) )
434  {
435  if ( box2df_left(a, b) )
436  return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
437  if ( box2df_right(a, b) )
438  return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
439  else
440  return (double)a->ymin - (double)b->ymax;
441  }
442  if ( box2df_below(a, b) )
443  {
444  if ( box2df_left(a, b) )
445  return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
446  if ( box2df_right(a, b) )
447  return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
448  else
449  return (double)b->ymin - (double)a->ymax;
450  }
451 
452  return FLT_MAX;
453 }
454 
455 static inline uint64_t
456 box2df_get_sortable_hash(const BOX2DF *b)
457 {
458  union floatuint {
459  uint32_t u;
460  float f;
461  } x, y;
462 
463  x.f = (b->xmax + b->xmin) / 2;
464  y.f = (b->ymax + b->ymin) / 2;
465 
466  return uint32_hilbert(y.u, x.u);
467 }
468 
474 int
475 gserialized_datum_get_internals_p(Datum gsdatum, GBOX *gbox, lwflags_t *flags, uint8_t *type, int32_t *srid)
476 {
477  int result = LW_SUCCESS;
478  GSERIALIZED *gpart = NULL;
479  int need_detoast = PG_GSERIALIZED_DATUM_NEEDS_DETOAST((struct varlena *)gsdatum);
480  if (need_detoast)
481  {
482  gpart = (GSERIALIZED *)PG_DETOAST_DATUM_SLICE(gsdatum, 0, gserialized_max_header_size());
483  }
484  else
485  {
486  gpart = (GSERIALIZED *)gsdatum;
487  }
488 
489  if (!gserialized_has_bbox(gpart) && need_detoast && LWSIZE_GET(gpart->size) >= gserialized_max_header_size())
490  {
491  /* The headers don't contain a bbox and the object is larger than what we retrieved, so
492  * we now detoast it completely */
493  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
494  gpart = (GSERIALIZED *)PG_DETOAST_DATUM(gsdatum);
495  }
496 
497  result = gserialized_get_gbox_p(gpart, gbox);
498  *flags = gserialized_get_lwflags(gpart);
499  *srid = gserialized_get_srid(gpart);
500  *type = gserialized_get_type(gpart);
501 
502  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
503  return result;
504 }
505 
513 int
514 gserialized_datum_get_gbox_p(Datum gsdatum, GBOX *gbox)
515 {
516  uint8_t type;
517  int32_t srid;
518  lwflags_t flags;
519 
520  return gserialized_datum_get_internals_p(gsdatum, gbox, &flags, &type, &srid);
521 }
522 
523 /* Note that this duplicates code from gserialized_datum_get_internals_p. It does so because
524  * accessing only the BOX2DF (as floats) is faster and this is a hot path function in indexes
525  */
526 int
527 gserialized_datum_get_box2df_p(Datum gsdatum, BOX2DF *box2df)
528 {
529  int result = LW_SUCCESS;
530  GSERIALIZED *gpart = NULL;
531  int need_detoast = PG_GSERIALIZED_DATUM_NEEDS_DETOAST((struct varlena *)gsdatum);
532  if (need_detoast)
533  {
534  gpart = (GSERIALIZED *)PG_DETOAST_DATUM_SLICE(gsdatum, 0, gserialized_max_header_size());
535  }
536  else
537  {
538  gpart = (GSERIALIZED *)gsdatum;
539  }
540 
541  if (gserialized_has_bbox(gpart))
542  {
543  size_t box_ndims;
544  const float *f = gserialized_get_float_box_p(gpart, &box_ndims);
545  memcpy(box2df, f, sizeof(BOX2DF));
546  result = LW_SUCCESS;
547  }
548  else
549  {
550  GBOX gbox = {0};
551  if (need_detoast && LWSIZE_GET(gpart->size) >= gserialized_max_header_size())
552  {
553  /* The headers don't contain a bbox and the object is larger than what we retrieved, so
554  * we now detoast it completely and recheck */
555  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
556  gpart = (GSERIALIZED *)PG_DETOAST_DATUM(gsdatum);
557  }
558  result = gserialized_get_gbox_p(gpart, &gbox);
559  if (result == LW_SUCCESS)
560  {
561  result = box2df_from_gbox_p(&gbox, box2df);
562  }
563  else {
564  /* Oh dear, has someone fed us an empty? */
565  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
566  gpart = (GSERIALIZED *)PG_DETOAST_DATUM(gsdatum);
567  if (gserialized_is_empty(gpart))
568  {
569  box2df_set_empty(box2df);
570  return LW_SUCCESS;
571  }
572  }
573  }
574 
575  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
576  return result;
577 }
578 
583 static int
584 gserialized_datum_predicate_2d(Datum gs1, Datum gs2, box2df_predicate predicate)
585 {
586  BOX2DF b1, b2, *br1=NULL, *br2=NULL;
587  POSTGIS_DEBUG(3, "entered function");
588 
589  if (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) br1 = &b1;
590  if (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) br2 = &b2;
591 
592  if ( predicate(br1, br2) )
593  {
594  POSTGIS_DEBUGF(3, "got boxes %s and %s", br1 ? box2df_to_string(&b1) : "(null)", br2 ? box2df_to_string(&b2) : "(null)");
595  return LW_TRUE;
596  }
597  return LW_FALSE;
598 }
599 
600 static int
601 gserialized_datum_predicate_box2df_geom_2d(const BOX2DF *br1, Datum gs2, box2df_predicate predicate)
602 {
603  BOX2DF b2, *br2=NULL;
604  POSTGIS_DEBUG(3, "entered function");
605 
606  if (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) br2 = &b2;
607 
608  if ( predicate(br1, br2) )
609  {
610  POSTGIS_DEBUGF(3, "got boxes %s", br2 ? box2df_to_string(&b2) : "(null)");
611  return LW_TRUE;
612  }
613  return LW_FALSE;
614 }
615 
616 /***********************************************************************
617 * BRIN 2-D Index Operator Functions
618 */
619 
621 Datum gserialized_contains_box2df_geom_2d(PG_FUNCTION_ARGS)
622 {
623  POSTGIS_DEBUG(3, "entered function");
624  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_contains) == LW_TRUE )
625  PG_RETURN_BOOL(true);
626 
627  PG_RETURN_BOOL(false);
628 }
629 
632 {
633  if ( box2df_contains((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
634  PG_RETURN_BOOL(true);
635 
636  PG_RETURN_BOOL(false);
637 }
638 
640 Datum gserialized_within_box2df_geom_2d(PG_FUNCTION_ARGS)
641 {
642  POSTGIS_DEBUG(3, "entered function");
643  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_within) == LW_TRUE )
644  PG_RETURN_BOOL(true);
645 
646  PG_RETURN_BOOL(false);
647 }
648 
650 Datum gserialized_within_box2df_box2df_2d(PG_FUNCTION_ARGS)
651 {
652  if ( box2df_within((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
653  PG_RETURN_BOOL(true);
654 
655  PG_RETURN_BOOL(false);
656 }
657 
659 Datum gserialized_overlaps_box2df_geom_2d(PG_FUNCTION_ARGS)
660 {
661  POSTGIS_DEBUG(3, "entered function");
662  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_overlaps) == LW_TRUE )
663  PG_RETURN_BOOL(true);
664 
665  PG_RETURN_BOOL(false);
666 }
667 
670 {
671  if ( box2df_overlaps((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
672  PG_RETURN_BOOL(true);
673 
674  PG_RETURN_BOOL(false);
675 }
676 
677 /***********************************************************************
678 * GiST 2-D Index Operator Functions
679 */
680 
682 Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS)
683 {
684  BOX2DF b1, b2;
685  Datum gs1 = PG_GETARG_DATUM(0);
686  Datum gs2 = PG_GETARG_DATUM(1);
687 
688  POSTGIS_DEBUG(3, "entered function");
689 
690  /* Must be able to build box for each argument (ie, not empty geometry). */
691  if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
693  {
694  double distance = box2df_distance_leaf_centroid(&b1, &b2);
695  POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
696  PG_RETURN_FLOAT8(distance);
697  }
698  PG_RETURN_FLOAT8(FLT_MAX);
699 }
700 
702 Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS)
703 {
704  BOX2DF b1, b2;
705  Datum gs1 = PG_GETARG_DATUM(0);
706  Datum gs2 = PG_GETARG_DATUM(1);
707 
708  POSTGIS_DEBUG(3, "entered function");
709 
710  /* Must be able to build box for each argument (ie, not empty geometry). */
711  if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
713  {
714  double distance = box2df_distance(&b1, &b2);
715  POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
716  PG_RETURN_FLOAT8(distance);
717  }
718  PG_RETURN_FLOAT8(FLT_MAX);
719 }
720 
722 Datum gserialized_same_2d(PG_FUNCTION_ARGS)
723 {
724  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_equals) == LW_TRUE )
725  PG_RETURN_BOOL(true);
726 
727  PG_RETURN_BOOL(false);
728 }
729 
731 Datum gserialized_within_2d(PG_FUNCTION_ARGS)
732 {
733  POSTGIS_DEBUG(3, "entered function");
734  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_within) == LW_TRUE )
735  PG_RETURN_BOOL(true);
736 
737  PG_RETURN_BOOL(false);
738 }
739 
741 Datum gserialized_contains_2d(PG_FUNCTION_ARGS)
742 {
743  POSTGIS_DEBUG(3, "entered function");
744  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_contains) == LW_TRUE )
745  PG_RETURN_BOOL(true);
746 
747  PG_RETURN_BOOL(false);
748 }
749 
751 Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS)
752 {
753  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overlaps) == LW_TRUE )
754  PG_RETURN_BOOL(true);
755 
756  PG_RETURN_BOOL(false);
757 }
758 
760 Datum gserialized_left_2d(PG_FUNCTION_ARGS)
761 {
762  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_left) == LW_TRUE )
763  PG_RETURN_BOOL(true);
764 
765  PG_RETURN_BOOL(false);
766 }
767 
769 Datum gserialized_right_2d(PG_FUNCTION_ARGS)
770 {
771  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_right) == LW_TRUE )
772  PG_RETURN_BOOL(true);
773 
774  PG_RETURN_BOOL(false);
775 }
776 
778 Datum gserialized_above_2d(PG_FUNCTION_ARGS)
779 {
780  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_above) == LW_TRUE )
781  PG_RETURN_BOOL(true);
782 
783  PG_RETURN_BOOL(false);
784 }
785 
787 Datum gserialized_below_2d(PG_FUNCTION_ARGS)
788 {
789  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_below) == LW_TRUE )
790  PG_RETURN_BOOL(true);
791 
792  PG_RETURN_BOOL(false);
793 }
794 
796 Datum gserialized_overleft_2d(PG_FUNCTION_ARGS)
797 {
798  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overleft) == LW_TRUE )
799  PG_RETURN_BOOL(true);
800 
801  PG_RETURN_BOOL(false);
802 }
803 
805 Datum gserialized_overright_2d(PG_FUNCTION_ARGS)
806 {
807  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overright) == LW_TRUE )
808  PG_RETURN_BOOL(true);
809 
810  PG_RETURN_BOOL(false);
811 }
812 
814 Datum gserialized_overabove_2d(PG_FUNCTION_ARGS)
815 {
816  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overabove) == LW_TRUE )
817  PG_RETURN_BOOL(true);
818 
819  PG_RETURN_BOOL(false);
820 }
821 
823 Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS)
824 {
825  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overbelow) == LW_TRUE )
826  PG_RETURN_BOOL(true);
827 
828  PG_RETURN_BOOL(false);
829 }
830 
831 
832 /***********************************************************************
833 * GiST Index Support Functions
834 */
835 
836 /*
837 ** GiST support function. Given a geography, return a "compressed"
838 ** version. In this case, we convert the geography into a geocentric
839 ** bounding box. If the geography already has the box embedded in it
840 ** we pull that out and hand it back.
841 */
843 Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS)
844 {
845  GISTENTRY *entry_in = (GISTENTRY*)PG_GETARG_POINTER(0);
846  GISTENTRY *entry_out = NULL;
847  BOX2DF bbox_out;
848  int result = LW_SUCCESS;
849 
850  POSTGIS_DEBUG(4, "[GIST] 'compress' function called");
851 
852  /*
853  ** Not a leaf key? There's nothing to do.
854  ** Return the input unchanged.
855  */
856  if ( ! entry_in->leafkey )
857  {
858  POSTGIS_DEBUG(4, "[GIST] non-leafkey entry, returning input unaltered");
859  PG_RETURN_POINTER(entry_in);
860  }
861 
862  POSTGIS_DEBUG(4, "[GIST] processing leafkey input");
863  entry_out = palloc(sizeof(GISTENTRY));
864 
865  /*
866  ** Null key? Make a copy of the input entry and
867  ** return.
868  */
869  if ( DatumGetPointer(entry_in->key) == NULL )
870  {
871  POSTGIS_DEBUG(4, "[GIST] leafkey is null");
872  gistentryinit(*entry_out, (Datum) 0, entry_in->rel,
873  entry_in->page, entry_in->offset, false);
874  POSTGIS_DEBUG(4, "[GIST] returning copy of input");
875  PG_RETURN_POINTER(entry_out);
876  }
877 
878  /* Extract our index key from the GiST entry. */
879  result = gserialized_datum_get_box2df_p(entry_in->key, &bbox_out);
880 
881  if (result == LW_FAILURE)
882  {
883  /* Treat failure like NULL */
884  POSTGIS_DEBUG(4, "[GIST] leafkey is null");
885  gistentryinit(*entry_out, (Datum) 0, entry_in->rel,
886  entry_in->page, entry_in->offset, false);
887  POSTGIS_DEBUG(4, "[GIST] returning copy of input");
888  PG_RETURN_POINTER(entry_out);
889  }
890 
891  /* Is the bounding box empty? Done! */
892  if (box2df_is_empty(&bbox_out))
893  {
894  gistentryinit(*entry_out, PointerGetDatum(box2df_copy(&bbox_out)),
895  entry_in->rel, entry_in->page, entry_in->offset, false);
896 
897  POSTGIS_DEBUG(4, "[GIST] empty geometry!");
898  PG_RETURN_POINTER(entry_out);
899  }
900 
901  POSTGIS_DEBUGF(4, "[GIST] got entry_in->key: %s", box2df_to_string(&bbox_out));
902 
903  /* Check all the dimensions for finite values */
904  if ( ! isfinite(bbox_out.xmax) || ! isfinite(bbox_out.xmin) ||
905  ! isfinite(bbox_out.ymax) || ! isfinite(bbox_out.ymin) )
906  {
907  box2df_set_finite(&bbox_out);
908  gistentryinit(*entry_out, PointerGetDatum(box2df_copy(&bbox_out)),
909  entry_in->rel, entry_in->page, entry_in->offset, false);
910 
911  POSTGIS_DEBUG(4, "[GIST] infinite geometry!");
912  PG_RETURN_POINTER(entry_out);
913  }
914 
915  /* Enure bounding box has minimums below maximums. */
916  box2df_validate(&bbox_out);
917 
918  /* Prepare GISTENTRY for return. */
919  gistentryinit(*entry_out, PointerGetDatum(box2df_copy(&bbox_out)),
920  entry_in->rel, entry_in->page, entry_in->offset, false);
921 
922  /* Return GISTENTRY. */
923  POSTGIS_DEBUG(4, "[GIST] 'compress' function complete");
924  PG_RETURN_POINTER(entry_out);
925 }
926 
927 
928 /*
929 ** GiST support function.
930 ** Decompress an entry. Unused for geography, so we return.
931 */
933 Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS)
934 {
935  POSTGIS_DEBUG(5, "[GIST] 'decompress' function called");
936  /* We don't decompress. Just return the input. */
937  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
938 }
939 
940 
941 
942 /*
943 ** GiST support function. Called from gserialized_gist_consistent below.
944 */
945 static inline bool gserialized_gist_consistent_leaf_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
946 {
947  bool retval;
948 
949  POSTGIS_DEBUGF(4, "[GIST] leaf consistent, strategy [%d], count[%i]",
950  strategy, g2d_counter_leaf++);
951 
952  switch (strategy)
953  {
954 
955  /* Basic overlaps */
956  case RTOverlapStrategyNumber:
957  retval = (bool) box2df_overlaps(key, query);
958  break;
959  case RTSameStrategyNumber:
960  retval = (bool) box2df_equals(key, query);
961  break;
962  case RTContainsStrategyNumber:
963  case RTOldContainsStrategyNumber:
964  retval = (bool) box2df_contains(key, query);
965  break;
966  case RTContainedByStrategyNumber:
967  case RTOldContainedByStrategyNumber:
968  retval = (bool) box2df_within(key, query);
969  break;
970 
971  /* To one side */
972  case RTAboveStrategyNumber:
973  retval = (bool) box2df_above(key, query);
974  break;
975  case RTBelowStrategyNumber:
976  retval = (bool) box2df_below(key, query);
977  break;
978  case RTRightStrategyNumber:
979  retval = (bool) box2df_right(key, query);
980  break;
981  case RTLeftStrategyNumber:
982  retval = (bool) box2df_left(key, query);
983  break;
984 
985  /* Overlapping to one side */
986  case RTOverAboveStrategyNumber:
987  retval = (bool) box2df_overabove(key, query);
988  break;
989  case RTOverBelowStrategyNumber:
990  retval = (bool) box2df_overbelow(key, query);
991  break;
992  case RTOverRightStrategyNumber:
993  retval = (bool) box2df_overright(key, query);
994  break;
995  case RTOverLeftStrategyNumber:
996  retval = (bool) box2df_overleft(key, query);
997  break;
998 
999  default:
1000  retval = false;
1001  }
1002 
1003  return (retval);
1004 }
1005 
1006 /*
1007 ** GiST support function. Called from gserialized_gist_consistent below.
1008 */
1009 static inline bool gserialized_gist_consistent_internal_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
1010 {
1011  bool retval;
1012 
1013  POSTGIS_DEBUGF(4, "[GIST] internal consistent, strategy [%d], count[%i], query[%s], key[%s]",
1014  strategy, g2d_counter_internal++, box2df_to_string(query), box2df_to_string(key) );
1015 
1016  switch (strategy)
1017  {
1018 
1019  /* Basic overlaps */
1020  case RTOverlapStrategyNumber:
1021  retval = (bool) box2df_overlaps(key, query);
1022  break;
1023  case RTSameStrategyNumber:
1024  case RTContainsStrategyNumber:
1025  case RTOldContainsStrategyNumber:
1026  retval = (bool) box2df_contains(key, query);
1027  break;
1028  case RTContainedByStrategyNumber:
1029  case RTOldContainedByStrategyNumber:
1030  retval = (bool) box2df_overlaps(key, query);
1031  break;
1032 
1033  /* To one side */
1034  case RTAboveStrategyNumber:
1035  retval = (bool)(!box2df_overbelow(key, query));
1036  break;
1037  case RTBelowStrategyNumber:
1038  retval = (bool)(!box2df_overabove(key, query));
1039  break;
1040  case RTRightStrategyNumber:
1041  retval = (bool)(!box2df_overleft(key, query));
1042  break;
1043  case RTLeftStrategyNumber:
1044  retval = (bool)(!box2df_overright(key, query));
1045  break;
1046 
1047  /* Overlapping to one side */
1048  case RTOverAboveStrategyNumber:
1049  retval = (bool)(!box2df_below(key, query));
1050  break;
1051  case RTOverBelowStrategyNumber:
1052  retval = (bool)(!box2df_above(key, query));
1053  break;
1054  case RTOverRightStrategyNumber:
1055  retval = (bool)(!box2df_left(key, query));
1056  break;
1057  case RTOverLeftStrategyNumber:
1058  retval = (bool)(!box2df_right(key, query));
1059  break;
1060 
1061  default:
1062  retval = false;
1063  }
1064 
1065  return (retval);
1066 }
1067 
1068 /*
1069 ** GiST support function. Take in a query and an entry and see what the
1070 ** relationship is, based on the query strategy.
1071 */
1073 Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS)
1074 {
1075  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1076  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1077  bool result;
1078  BOX2DF query_gbox_index;
1079 
1080  /* PostgreSQL 8.4 and later require the RECHECK flag to be set here,
1081  rather than being supplied as part of the operator class definition */
1082  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1083 
1084  /* We set recheck to false to avoid repeatedly pulling every "possibly matched" geometry
1085  out during index scans. For cases when the geometries are large, rechecking
1086  can make things twice as slow. */
1087  *recheck = false;
1088 
1089  POSTGIS_DEBUG(4, "[GIST] 'consistent' function called");
1090 
1091  /* Quick sanity check on query argument. */
1092  if ( DatumGetPointer(PG_GETARG_DATUM(1)) == NULL )
1093  {
1094  POSTGIS_DEBUG(4, "[GIST] null query pointer (!?!), returning false");
1095  PG_RETURN_BOOL(false); /* NULL query! This is screwy! */
1096  }
1097 
1098  /* Quick sanity check on entry key. */
1099  if ( DatumGetPointer(entry->key) == NULL )
1100  {
1101  POSTGIS_DEBUG(4, "[GIST] null index entry, returning false");
1102  PG_RETURN_BOOL(false); /* NULL entry! */
1103  }
1104 
1105  /* Null box should never make this far. */
1106  if ( gserialized_datum_get_box2df_p(PG_GETARG_DATUM(1), &query_gbox_index) == LW_FAILURE )
1107  {
1108  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1109  PG_RETURN_BOOL(false);
1110  }
1111 
1112  /* Treat leaf node tests different from internal nodes */
1113  if (GIST_LEAF(entry))
1114  {
1116  (BOX2DF*)DatumGetPointer(entry->key),
1117  &query_gbox_index, strategy);
1118  }
1119  else
1120  {
1122  (BOX2DF*)DatumGetPointer(entry->key),
1123  &query_gbox_index, strategy);
1124  }
1125 
1126  PG_RETURN_BOOL(result);
1127 }
1128 
1129 
1130 /*
1131 ** GiST support function. Take in a query and an entry and return the "distance"
1132 ** between them.
1133 **
1134 ** Given an index entry p and a query value q, this function determines the
1135 ** index entry's "distance" from the query value. This function must be
1136 ** supplied if the operator class contains any ordering operators. A query
1137 ** using the ordering operator will be implemented by returning index entries
1138 ** with the smallest "distance" values first, so the results must be consistent
1139 ** with the operator's semantics. For a leaf index entry the result just
1140 ** represents the distance to the index entry; for an internal tree node, the
1141 ** result must be the smallest distance that any child entry could have.
1142 **
1143 ** Strategy 13 = true distance tests <->
1144 ** Strategy 14 = box-based distance tests <#>
1145 */
1147 Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
1148 {
1149  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1150  BOX2DF query_box;
1151  BOX2DF *entry_box;
1152  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1153  double distance;
1154  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1155 
1156  POSTGIS_DEBUG(4, "[GIST] 'distance' function called");
1157 
1158  /* We are using '13' as the gist true-distance <-> strategy number
1159  * and '14' as the gist distance-between-boxes <#> strategy number */
1160  if ( strategy != 13 && strategy != 14 ) {
1161  elog(ERROR, "unrecognized strategy number: %d", strategy);
1162  PG_RETURN_FLOAT8(FLT_MAX);
1163  }
1164 
1165  /* Null box should never make this far. */
1166  if ( gserialized_datum_get_box2df_p(PG_GETARG_DATUM(1), &query_box) == LW_FAILURE )
1167  {
1168  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1169  PG_RETURN_FLOAT8(FLT_MAX);
1170  }
1171 
1172  /* Get the entry box */
1173  entry_box = (BOX2DF*)DatumGetPointer(entry->key);
1174 
1175  /* Box-style distance test */
1176  if ( strategy == 14 ) /* operator <#> */
1177  {
1178  distance = box2df_distance(entry_box, &query_box);
1179  }
1180  /* True distance test (formerly centroid distance) */
1181  else if ( strategy == 13 ) /* operator <-> */
1182  {
1183  /* In all cases, since we only have keys (boxes) we'll return */
1184  /* the minimum possible distance, which is the box2df_distance */
1185  /* and let the recheck sort things out in the case of leaves */
1186  distance = box2df_distance(entry_box, &query_box);
1187 
1188  if (GIST_LEAF(entry))
1189  *recheck = true;
1190  }
1191  else
1192  {
1193  elog(ERROR, "%s: reached unreachable code", __func__);
1194  PG_RETURN_NULL();
1195  }
1196 
1197  PG_RETURN_FLOAT8(distance);
1198 }
1199 
1200 /*
1201 ** Function to pack floats of different realms.
1202 ** This function serves to pack bit flags inside float type.
1203 ** Result value represent can be from two different "realms".
1204 ** Every value from realm 1 is greater than any value from realm 0.
1205 ** Values from the same realm loose one bit of precision.
1206 **
1207 ** This technique is possible due to floating point numbers specification
1208 ** according to IEEE 754: exponent bits are highest
1209 ** (excluding sign bits, but here penalty is always positive). If float a is
1210 ** greater than float b, integer A with same bit representation as a is greater
1211 ** than integer B with same bits as b.
1212 */
1213 static inline float
1214 pack_float(const float value, const uint8_t realm)
1215 {
1216  union {
1217  float f;
1218  struct {
1219  unsigned value : 31, sign : 1;
1220  } vbits;
1221  struct {
1222  unsigned value : 30, realm : 1, sign : 1;
1223  } rbits;
1224  } a;
1225 
1226  a.f = value;
1227  a.rbits.value = a.vbits.value >> 1;
1228  a.rbits.realm = realm;
1229 
1230  return a.f;
1231 }
1232 
1233 static inline float
1234 box2df_penalty(const BOX2DF *b1, const BOX2DF *b2)
1235 {
1236  float b1xmin = b1->xmin, b1xmax = b1->xmax;
1237  float b1ymin = b1->ymin, b1ymax = b1->ymax;
1238  float b2xmin = b2->xmin, b2xmax = b2->xmax;
1239  float b2ymin = b2->ymin, b2ymax = b2->ymax;
1240 
1241  float box_union_xmin = Min(b1xmin, b2xmin), box_union_xmax = Max(b1xmax, b2xmax);
1242  float box_union_ymin = Min(b1ymin, b2ymin), box_union_ymax = Max(b1ymax, b2ymax);
1243 
1244  float b1dx = b1xmax - b1xmin, b1dy = b1ymax - b1ymin;
1245  float box_union_dx = box_union_xmax - box_union_xmin, box_union_dy = box_union_ymax - box_union_ymin;
1246 
1247  float box_union_area = box_union_dx * box_union_dy, box1area = b1dx * b1dy;
1248  float box_union_edge = box_union_dx + box_union_dy, box1edge = b1dx + b1dy;
1249 
1250  float area_extension = box_union_area - box1area;
1251  float edge_extension = box_union_edge - box1edge;
1252 
1253  /* REALM 1: Area extension is nonzero, return it */
1254  if (area_extension > FLT_EPSILON)
1255  return pack_float(area_extension, 1);
1256  /* REALM 0: Area extension is zero, return nonzero edge extension */
1257  else if (edge_extension > FLT_EPSILON)
1258  return pack_float(edge_extension, 0);
1259 
1260  return 0;
1261 }
1262 
1263 static inline float
1264 box2df_penalty_single(const BOX2DF *box)
1265 {
1266  float boxxmin = box->xmin, boxxmax = box->xmax;
1267  float boxymin = box->ymin, boxymax = box->ymax;
1268 
1269  float dx = boxxmax - boxxmin, dy = boxymax - boxymin;
1270 
1271  float area = dx * dy;
1272  float edge = dx + dy;
1273 
1274  /* REALM 1: Area is nonzero, return it */
1275  if (area > FLT_EPSILON)
1276  return pack_float(area, 1);
1277  /* REALM 0: Area is zero, return edge */
1278  else if (edge > FLT_EPSILON)
1279  return pack_float(edge, 0);
1280 
1281  return 0;
1282 }
1283 
1284 /*
1285 ** GiST support function. Calculate the "penalty" cost of adding this entry into an existing entry.
1286 ** Calculate the change in volume of the old entry once the new entry is added.
1287 */
1289 Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS)
1290 {
1291  GISTENTRY *origentry = (GISTENTRY*) PG_GETARG_POINTER(0);
1292  GISTENTRY *newentry = (GISTENTRY*) PG_GETARG_POINTER(1);
1293  float *result = (float*) PG_GETARG_POINTER(2);
1294  BOX2DF *b1, *b2;
1295 
1296  b1 = (BOX2DF *)DatumGetPointer(origentry->key);
1297  b2 = (BOX2DF *)DatumGetPointer(newentry->key);
1298 
1299  /* Penalty value of 0 has special code path in Postgres's gistchoose.
1300  * It is used as an early exit condition in subtree loop, allowing faster tree descend.
1301  * For multi-column index, it lets next column break the tie, possibly more confidently.
1302  */
1303  *result = 0;
1304 
1305  if (b1 && b2 && !box2df_is_empty(b1) && !box2df_is_empty(b2))
1306  *result = box2df_penalty(b1, b2);
1307 
1308  PG_RETURN_POINTER(result);
1309 }
1310 
1311 /*
1312 ** GiST support function. Merge all the boxes in a page into one master box.
1313 */
1315 Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS)
1316 {
1317  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
1318  int *sizep = (int *) PG_GETARG_POINTER(1); /* Size of the return value */
1319  int numranges, i;
1320  BOX2DF *box_cur, *box_union;
1321 
1322  POSTGIS_DEBUG(4, "[GIST] 'union' function called");
1323 
1324  numranges = entryvec->n;
1325 
1326  box_cur = (BOX2DF*) DatumGetPointer(entryvec->vector[0].key);
1327 
1328  box_union = box2df_copy(box_cur);
1329 
1330  for ( i = 1; i < numranges; i++ )
1331  {
1332  box_cur = (BOX2DF*) DatumGetPointer(entryvec->vector[i].key);
1333  box2df_merge(box_union, box_cur);
1334  }
1335 
1336  *sizep = sizeof(BOX2DF);
1337 
1338  POSTGIS_DEBUGF(4, "[GIST] 'union', numranges(%i), pageunion %s", numranges, box2df_to_string(box_union));
1339 
1340  PG_RETURN_POINTER(box_union);
1341 
1342 }
1343 
1344 /*
1345 ** GiST support function. Test equality of keys.
1346 */
1348 Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS)
1349 {
1350  BOX2DF *b1 = (BOX2DF*)PG_GETARG_POINTER(0);
1351  BOX2DF *b2 = (BOX2DF*)PG_GETARG_POINTER(1);
1352  bool *result = (bool*)PG_GETARG_POINTER(2);
1353 
1354  POSTGIS_DEBUG(4, "[GIST] 'same' function called");
1355 
1356  *result = box2df_equals(b1, b2);
1357 
1358  PG_RETURN_POINTER(result);
1359 }
1360 
1361 static int
1362 gserialized_gist_cmp_abbrev_2d(Datum x, Datum y, SortSupport ssup)
1363 {
1364  /* Empty is a special case */
1365  if (x == 0 || y == 0 || x == y)
1366  return 0; /* 0 means "ask bigger comparator" and not equality*/
1367  else if (x > y)
1368  return 1;
1369  else
1370  return -1;
1371 }
1372 
1373 static bool
1374 gserialized_gist_abbrev_abort_2d(int memtupcount, SortSupport ssup)
1375 {
1376  return LW_FALSE;
1377 }
1378 
1379 static Datum
1380 gserialized_gist_abbrev_convert_2d(Datum original, SortSupport ssup)
1381 {
1382  return box2df_get_sortable_hash((BOX2DF *)original);
1383 }
1384 
1385 static int
1386 gserialized_gist_cmp_full_2d(Datum a, Datum b, SortSupport ssup)
1387 {
1388  BOX2DF *b1 = (BOX2DF *)a;
1389  BOX2DF *b2 = (BOX2DF *)b;
1390  uint64_t hash1, hash2;
1391  int cmp;
1392 
1393  cmp = memcmp(b1, b2, sizeof(BOX2DF));
1394  if (cmp == 0)
1395  return 0;
1396 
1397  hash1 = box2df_get_sortable_hash(b1);
1398  hash2 = box2df_get_sortable_hash(b2);
1399  if (hash1 > hash2)
1400  return 1;
1401  else if (hash1 < hash2)
1402  return -1;
1403 
1404  return cmp > 0 ? 1 : -1;
1405 }
1406 
1408 Datum gserialized_gist_sortsupport_2d(PG_FUNCTION_ARGS)
1409 {
1410  SortSupport ssup = (SortSupport)PG_GETARG_POINTER(0);
1411 
1412  ssup->comparator = gserialized_gist_cmp_full_2d;
1413  ssup->ssup_extra = NULL;
1414  /* Enable sortsupport only on 64 bit Datum */
1415  if (ssup->abbreviate && sizeof(Datum) == 8)
1416  {
1417  ssup->comparator = gserialized_gist_cmp_abbrev_2d;
1418  ssup->abbrev_converter = gserialized_gist_abbrev_convert_2d;
1419  ssup->abbrev_abort = gserialized_gist_abbrev_abort_2d;
1420  ssup->abbrev_full_comparator = gserialized_gist_cmp_full_2d;
1421  }
1422 
1423  PG_RETURN_VOID();
1424 }
1425 
1426 /*
1427  * Adjust BOX2DF b boundaries with insertion of addon.
1428  */
1429 static void
1430 adjustBox(BOX2DF *b, BOX2DF *addon)
1431 {
1432  if (b->xmax < addon->xmax || isnan(b->xmax))
1433  b->xmax = addon->xmax;
1434  if (b->xmin > addon->xmin || isnan(b->xmin))
1435  b->xmin = addon->xmin;
1436  if (b->ymax < addon->ymax || isnan(b->ymax))
1437  b->ymax = addon->ymax;
1438  if (b->ymin > addon->ymin || isnan(b->ymin))
1439  b->ymin = addon->ymin;
1440 }
1441 
1442 /*
1443  * Trivial split: half of entries will be placed on one page
1444  * and another half - to another
1445  */
1446 static void
1447 fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
1448 {
1449  OffsetNumber i,
1450  maxoff;
1451  BOX2DF *unionL = NULL,
1452  *unionR = NULL;
1453  int nbytes;
1454 
1455  maxoff = entryvec->n - 1;
1456 
1457  nbytes = (maxoff + 2) * sizeof(OffsetNumber);
1458  v->spl_left = (OffsetNumber *) palloc(nbytes);
1459  v->spl_right = (OffsetNumber *) palloc(nbytes);
1460  v->spl_nleft = v->spl_nright = 0;
1461 
1462  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1463  {
1464  BOX2DF *cur = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
1465 
1466  if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
1467  {
1468  v->spl_left[v->spl_nleft] = i;
1469  if (unionL == NULL)
1470  {
1471  unionL = (BOX2DF *) palloc(sizeof(BOX2DF));
1472  *unionL = *cur;
1473  }
1474  else
1475  adjustBox(unionL, cur);
1476 
1477  v->spl_nleft++;
1478  }
1479  else
1480  {
1481  v->spl_right[v->spl_nright] = i;
1482  if (unionR == NULL)
1483  {
1484  unionR = (BOX2DF *) palloc(sizeof(BOX2DF));
1485  *unionR = *cur;
1486  }
1487  else
1488  adjustBox(unionR, cur);
1489 
1490  v->spl_nright++;
1491  }
1492  }
1493 
1494  if (v->spl_ldatum_exists)
1495  adjustBox(unionL, (BOX2DF *) DatumGetPointer(v->spl_ldatum));
1496  v->spl_ldatum = PointerGetDatum(unionL);
1497 
1498  if (v->spl_rdatum_exists)
1499  adjustBox(unionR, (BOX2DF *) DatumGetPointer(v->spl_rdatum));
1500  v->spl_rdatum = PointerGetDatum(unionR);
1501 
1502  v->spl_ldatum_exists = v->spl_rdatum_exists = false;
1503 }
1504 
1505 /*
1506  * Represents information about an entry that can be placed to either group
1507  * without affecting overlap over selected axis ("common entry").
1508  */
1509 typedef struct
1510 {
1511  /* Index of entry in the initial array */
1512  int index;
1513  /* Delta between penalties of entry insertion into different groups */
1514  float delta;
1515 } CommonEntry;
1516 
1517 /*
1518  * Context for g_box_consider_split. Contains information about currently
1519  * selected split and some general information.
1520  */
1521 typedef struct
1522 {
1523  int entriesCount; /* total number of entries being split */
1524  BOX2DF boundingBox; /* minimum bounding box across all entries */
1525 
1526  /* Information about currently selected split follows */
1527 
1528  bool first; /* true if no split was selected yet */
1529 
1530  float leftUpper; /* upper bound of left interval */
1531  float rightLower; /* lower bound of right interval */
1532 
1533  float4 ratio;
1534  float4 overlap;
1535  int dim; /* axis of this split */
1536  float range; /* width of general MBR projection to the
1537  * selected axis */
1539 
1540 /*
1541  * Interval represents projection of box to axis.
1542  */
1543 typedef struct
1544 {
1545  float lower,
1547 } SplitInterval;
1548 
1549 /*
1550  * Interval comparison function by lower bound of the interval;
1551  */
1552 static int
1553 interval_cmp_lower(const void *i1, const void *i2)
1554 {
1555  float lower1 = ((const SplitInterval *) i1)->lower,
1556  lower2 = ((const SplitInterval *) i2)->lower;
1557 
1558  if (isnan(lower1))
1559  {
1560  if (isnan(lower2))
1561  return 0;
1562  else
1563  return 1;
1564  }
1565  else if (isnan(lower2))
1566  {
1567  return -1;
1568  }
1569  else
1570  {
1571  if (lower1 < lower2)
1572  return -1;
1573  else if (lower1 > lower2)
1574  return 1;
1575  else
1576  return 0;
1577  }
1578 }
1579 
1580 
1581 /*
1582  * Interval comparison function by upper bound of the interval;
1583  */
1584 static int
1585 interval_cmp_upper(const void *i1, const void *i2)
1586 {
1587  float upper1 = ((const SplitInterval *) i1)->upper,
1588  upper2 = ((const SplitInterval *) i2)->upper;
1589 
1590  if (isnan(upper1))
1591  {
1592  if (isnan(upper2))
1593  return 0;
1594  else
1595  return -1;
1596  }
1597  else if (isnan(upper2))
1598  {
1599  return 1;
1600  }
1601  else
1602  {
1603  if (upper1 < upper2)
1604  return -1;
1605  else if (upper1 > upper2)
1606  return 1;
1607  else
1608  return 0;
1609  }
1610 }
1611 
1612 /*
1613  * Replace negative value with zero.
1614  */
1615 static inline float
1616 non_negative(float val)
1617 {
1618  if (val >= 0.0f)
1619  return val;
1620  else
1621  return 0.0f;
1622 }
1623 
1624 /*
1625  * Consider replacement of currently selected split with the better one.
1626  */
1627 static inline void
1629  float rightLower, int minLeftCount,
1630  float leftUpper, int maxLeftCount)
1631 {
1632  int leftCount,
1633  rightCount;
1634  float4 ratio,
1635  overlap;
1636  float range;
1637 
1638  POSTGIS_DEBUGF(5, "consider split: dimNum = %d, rightLower = %f, "
1639  "minLeftCount = %d, leftUpper = %f, maxLeftCount = %d ",
1640  dimNum, rightLower, minLeftCount, leftUpper, maxLeftCount);
1641 
1642  /*
1643  * Calculate entries distribution ratio assuming most uniform distribution
1644  * of common entries.
1645  */
1646  if (minLeftCount >= (context->entriesCount + 1) / 2)
1647  {
1648  leftCount = minLeftCount;
1649  }
1650  else
1651  {
1652  if (maxLeftCount <= context->entriesCount / 2)
1653  leftCount = maxLeftCount;
1654  else
1655  leftCount = context->entriesCount / 2;
1656  }
1657  rightCount = context->entriesCount - leftCount;
1658 
1659  /*
1660  * Ratio of split - quotient between size of lesser group and total
1661  * entries count.
1662  */
1663  ratio = ((float4) Min(leftCount, rightCount)) /
1664  ((float4) context->entriesCount);
1665 
1666  if (ratio > LIMIT_RATIO)
1667  {
1668  bool selectthis = false;
1669 
1670  /*
1671  * The ratio is acceptable, so compare current split with previously
1672  * selected one. Between splits of one dimension we search for minimal
1673  * overlap (allowing negative values) and minimal ration (between same
1674  * overlaps. We switch dimension if find less overlap (non-negative)
1675  * or less range with same overlap.
1676  */
1677  if (dimNum == 0)
1678  range = context->boundingBox.xmax - context->boundingBox.xmin;
1679  else
1680  range = context->boundingBox.ymax - context->boundingBox.ymin;
1681 
1682  overlap = (leftUpper - rightLower) / range;
1683 
1684  /* If there is no previous selection, select this */
1685  if (context->first)
1686  selectthis = true;
1687  else if (context->dim == dimNum)
1688  {
1689  /*
1690  * Within the same dimension, choose the new split if it has a
1691  * smaller overlap, or same overlap but better ratio.
1692  */
1693  if (overlap < context->overlap ||
1694  (overlap == context->overlap && ratio > context->ratio))
1695  selectthis = true;
1696  }
1697  else
1698  {
1699  /*
1700  * Across dimensions, choose the new split if it has a smaller
1701  * *non-negative* overlap, or same *non-negative* overlap but
1702  * bigger range. This condition differs from the one described in
1703  * the article. On the datasets where leaf MBRs don't overlap
1704  * themselves, non-overlapping splits (i.e. splits which have zero
1705  * *non-negative* overlap) are frequently possible. In this case
1706  * splits tends to be along one dimension, because most distant
1707  * non-overlapping splits (i.e. having lowest negative overlap)
1708  * appears to be in the same dimension as in the previous split.
1709  * Therefore MBRs appear to be very prolonged along another
1710  * dimension, which leads to bad search performance. Using range
1711  * as the second split criteria makes MBRs more quadratic. Using
1712  * *non-negative* overlap instead of overlap as the first split
1713  * criteria gives to range criteria a chance to matter, because
1714  * non-overlapping splits are equivalent in this criteria.
1715  */
1716  if (non_negative(overlap) < non_negative(context->overlap) ||
1717  (range > context->range &&
1718  non_negative(overlap) <= non_negative(context->overlap)))
1719  selectthis = true;
1720  }
1721 
1722  if (selectthis)
1723  {
1724  /* save information about selected split */
1725  context->first = false;
1726  context->ratio = ratio;
1727  context->range = range;
1728  context->overlap = overlap;
1729  context->rightLower = rightLower;
1730  context->leftUpper = leftUpper;
1731  context->dim = dimNum;
1732  POSTGIS_DEBUG(5, "split selected");
1733  }
1734  }
1735 }
1736 
1737 /*
1738  * Compare common entries by their deltas.
1739  */
1740 static int
1741 common_entry_cmp(const void *i1, const void *i2)
1742 {
1743  float delta1 = ((const CommonEntry *) i1)->delta,
1744  delta2 = ((const CommonEntry *) i2)->delta;
1745 
1746  if (delta1 < delta2)
1747  return -1;
1748  else if (delta1 > delta2)
1749  return 1;
1750  else
1751  return 0;
1752 }
1753 
1754 /*
1755  * --------------------------------------------------------------------------
1756  * Double sorting split algorithm. This is used for both boxes and points.
1757  *
1758  * The algorithm finds split of boxes by considering splits along each axis.
1759  * Each entry is first projected as an interval on the X-axis, and different
1760  * ways to split the intervals into two groups are considered, trying to
1761  * minimize the overlap of the groups. Then the same is repeated for the
1762  * Y-axis, and the overall best split is chosen. The quality of a split is
1763  * determined by overlap along that axis and some other criteria (see
1764  * g_box_consider_split).
1765  *
1766  * After that, all the entries are divided into three groups:
1767  *
1768  * 1) Entries which should be placed to the left group
1769  * 2) Entries which should be placed to the right group
1770  * 3) "Common entries" which can be placed to any of groups without affecting
1771  * of overlap along selected axis.
1772  *
1773  * The common entries are distributed by minimizing penalty.
1774  *
1775  * For details see:
1776  * "A new double sorting-based node splitting algorithm for R-tree", A. Korotkov
1777  * http://syrcose.ispras.ru/2011/files/SYRCoSE2011_Proceedings.pdf#page=36
1778  * --------------------------------------------------------------------------
1779  */
1781 Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS)
1782 {
1783  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
1784  GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
1785  OffsetNumber i,
1786  maxoff;
1787  ConsiderSplitContext context;
1788  BOX2DF *leftBox, *rightBox;
1789  int dim,
1790  commonEntriesCount;
1791  SplitInterval *intervalsLower,
1792  *intervalsUpper;
1793  CommonEntry *commonEntries;
1794  int nentries;
1795 
1796  POSTGIS_DEBUG(3, "[GIST] 'picksplit' entered");
1797 
1798  memset(&context, 0, sizeof(ConsiderSplitContext));
1799 
1800  maxoff = entryvec->n - 1;
1801  nentries = context.entriesCount = maxoff - FirstOffsetNumber + 1;
1802 
1803  /* Allocate arrays for intervals along axes */
1804  intervalsLower = (SplitInterval *) palloc(nentries * sizeof(SplitInterval));
1805  intervalsUpper = (SplitInterval *) palloc(nentries * sizeof(SplitInterval));
1806 
1807  /*
1808  * Calculate the overall minimum bounding box over all the entries.
1809  */
1810  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1811  {
1812  BOX2DF *box = (BOX2DF *)DatumGetPointer(entryvec->vector[i].key);
1813  if (i == FirstOffsetNumber)
1814  context.boundingBox = *box;
1815  else
1816  adjustBox(&context.boundingBox, box);
1817  }
1818 
1819  POSTGIS_DEBUGF(4, "boundingBox is %s", box2df_to_string(
1820  &context.boundingBox));
1821 
1822  /*
1823  * Iterate over axes for optimal split searching.
1824  */
1825  context.first = true; /* nothing selected yet */
1826  for (dim = 0; dim < 2; dim++)
1827  {
1828  float leftUpper,
1829  rightLower;
1830  int i1,
1831  i2;
1832 
1833  /* Project each entry as an interval on the selected axis. */
1834  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1835  {
1836  BOX2DF *box = (BOX2DF *)DatumGetPointer(entryvec->vector[i].key);
1837  if (dim == 0)
1838  {
1839  intervalsLower[i - FirstOffsetNumber].lower = box->xmin;
1840  intervalsLower[i - FirstOffsetNumber].upper = box->xmax;
1841  }
1842  else
1843  {
1844  intervalsLower[i - FirstOffsetNumber].lower = box->ymin;
1845  intervalsLower[i - FirstOffsetNumber].upper = box->ymax;
1846  }
1847  }
1848 
1849  /*
1850  * Make two arrays of intervals: one sorted by lower bound and another
1851  * sorted by upper bound.
1852  */
1853  memcpy(intervalsUpper, intervalsLower,
1854  sizeof(SplitInterval) * nentries);
1855  qsort(intervalsLower, nentries, sizeof(SplitInterval),
1857  qsort(intervalsUpper, nentries, sizeof(SplitInterval),
1859 
1860  /*----
1861  * The goal is to form a left and right interval, so that every entry
1862  * interval is contained by either left or right interval (or both).
1863  *
1864  * For example, with the intervals (0,1), (1,3), (2,3), (2,4):
1865  *
1866  * 0 1 2 3 4
1867  * +-+
1868  * +---+
1869  * +-+
1870  * +---+
1871  *
1872  * The left and right intervals are of the form (0,a) and (b,4).
1873  * We first consider splits where b is the lower bound of an entry.
1874  * We iterate through all entries, and for each b, calculate the
1875  * smallest possible a. Then we consider splits where a is the
1876  * upper bound of an entry, and for each a, calculate the greatest
1877  * possible b.
1878  *
1879  * In the above example, the first loop would consider splits:
1880  * b=0: (0,1)-(0,4)
1881  * b=1: (0,1)-(1,4)
1882  * b=2: (0,3)-(2,4)
1883  *
1884  * And the second loop:
1885  * a=1: (0,1)-(1,4)
1886  * a=3: (0,3)-(2,4)
1887  * a=4: (0,4)-(2,4)
1888  */
1889 
1890  /*
1891  * Iterate over lower bound of right group, finding smallest possible
1892  * upper bound of left group.
1893  */
1894  i1 = 0;
1895  i2 = 0;
1896  rightLower = intervalsLower[i1].lower;
1897  leftUpper = intervalsUpper[i2].lower;
1898  while (true)
1899  {
1900  /*
1901  * Find next lower bound of right group.
1902  */
1903  while (i1 < nentries && (rightLower == intervalsLower[i1].lower ||
1904  isnan(intervalsLower[i1].lower)))
1905  {
1906  leftUpper = Max(leftUpper, intervalsLower[i1].upper);
1907  i1++;
1908  }
1909  if (i1 >= nentries)
1910  break;
1911  rightLower = intervalsLower[i1].lower;
1912 
1913  /*
1914  * Find count of intervals which anyway should be placed to the
1915  * left group.
1916  */
1917  while (i2 < nentries && intervalsUpper[i2].upper <= leftUpper)
1918  i2++;
1919 
1920  /*
1921  * Consider found split.
1922  */
1923  g_box_consider_split(&context, dim, rightLower, i1, leftUpper, i2);
1924  }
1925 
1926  /*
1927  * Iterate over upper bound of left group finding greatest possible
1928  * lower bound of right group.
1929  */
1930  i1 = nentries - 1;
1931  i2 = nentries - 1;
1932  rightLower = intervalsLower[i1].upper;
1933  leftUpper = intervalsUpper[i2].upper;
1934  while (true)
1935  {
1936  /*
1937  * Find next upper bound of left group.
1938  */
1939  while (i2 >= 0 && (leftUpper == intervalsUpper[i2].upper ||
1940  isnan(intervalsUpper[i2].upper)))
1941  {
1942  rightLower = Min(rightLower, intervalsUpper[i2].lower);
1943  i2--;
1944  }
1945  if (i2 < 0)
1946  break;
1947  leftUpper = intervalsUpper[i2].upper;
1948 
1949  /*
1950  * Find count of intervals which anyway should be placed to the
1951  * right group.
1952  */
1953  while (i1 >= 0 && intervalsLower[i1].lower >= rightLower)
1954  i1--;
1955 
1956  /*
1957  * Consider found split.
1958  */
1959  g_box_consider_split(&context, dim,
1960  rightLower, i1 + 1, leftUpper, i2 + 1);
1961  }
1962  }
1963 
1964  /*
1965  * If we failed to find any acceptable splits, use trivial split.
1966  */
1967  if (context.first)
1968  {
1969  POSTGIS_DEBUG(4, "no acceptable splits, trivial split");
1970  fallbackSplit(entryvec, v);
1971  PG_RETURN_POINTER(v);
1972  }
1973 
1974  /*
1975  * Ok, we have now selected the split across one axis.
1976  *
1977  * While considering the splits, we already determined that there will be
1978  * enough entries in both groups to reach the desired ratio, but we did
1979  * not memorize which entries go to which group. So determine that now.
1980  */
1981 
1982  POSTGIS_DEBUGF(4, "split direction: %d", context.dim);
1983 
1984  /* Allocate vectors for results */
1985  v->spl_left = (OffsetNumber *) palloc(nentries * sizeof(OffsetNumber));
1986  v->spl_right = (OffsetNumber *) palloc(nentries * sizeof(OffsetNumber));
1987  v->spl_nleft = 0;
1988  v->spl_nright = 0;
1989 
1990  /* Allocate bounding boxes of left and right groups */
1991  leftBox = palloc0(sizeof(BOX2DF));
1992  rightBox = palloc0(sizeof(BOX2DF));
1993 
1994  /*
1995  * Allocate an array for "common entries" - entries which can be placed to
1996  * either group without affecting overlap along selected axis.
1997  */
1998  commonEntriesCount = 0;
1999  commonEntries = (CommonEntry *) palloc(nentries * sizeof(CommonEntry));
2000 
2001  /* Helper macros to place an entry in the left or right group */
2002 #define PLACE_LEFT(box, off) \
2003  do { \
2004  if (v->spl_nleft > 0) \
2005  adjustBox(leftBox, box); \
2006  else \
2007  *leftBox = *(box); \
2008  v->spl_left[v->spl_nleft++] = off; \
2009  } while(0)
2010 
2011 #define PLACE_RIGHT(box, off) \
2012  do { \
2013  if (v->spl_nright > 0) \
2014  adjustBox(rightBox, box); \
2015  else \
2016  *rightBox = *(box); \
2017  v->spl_right[v->spl_nright++] = off; \
2018  } while(0)
2019 
2020  /*
2021  * Distribute entries which can be distributed unambiguously, and collect
2022  * common entries.
2023  */
2024  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
2025  {
2026  float lower,
2027  upper;
2028 
2029  /*
2030  * Get upper and lower bounds along selected axis.
2031  */
2032  BOX2DF *box = (BOX2DF *)DatumGetPointer(entryvec->vector[i].key);
2033  if (context.dim == 0)
2034  {
2035  lower = box->xmin;
2036  upper = box->xmax;
2037  }
2038  else
2039  {
2040  lower = box->ymin;
2041  upper = box->ymax;
2042  }
2043 
2044  if (upper <= context.leftUpper || isnan(upper))
2045  {
2046  /* Fits to the left group */
2047  if (lower >= context.rightLower || isnan(lower))
2048  {
2049  /* Fits also to the right group, so "common entry" */
2050  commonEntries[commonEntriesCount++].index = i;
2051  }
2052  else
2053  {
2054  /* Doesn't fit to the right group, so join to the left group */
2055  PLACE_LEFT(box, i);
2056  }
2057  }
2058  else
2059  {
2060  /*
2061  * Each entry should fit on either left or right group. Since this
2062  * entry didn't fit on the left group, it better fit in the right
2063  * group.
2064  */
2065  Assert(lower >= context.rightLower);
2066 
2067  /* Doesn't fit to the left group, so join to the right group */
2068  PLACE_RIGHT(box, i);
2069  }
2070  }
2071 
2072  POSTGIS_DEBUGF(4, "leftBox is %s", box2df_to_string(leftBox));
2073  POSTGIS_DEBUGF(4, "rightBox is %s", box2df_to_string(rightBox));
2074 
2075  /*
2076  * Distribute "common entries", if any.
2077  */
2078  if (commonEntriesCount > 0)
2079  {
2080  /* Calculate minimum number of entries that must be placed in both groups, to reach LIMIT_RATIO. */
2081  int m = ceil(LIMIT_RATIO * (double)nentries);
2082 
2083  /* Recursive picksplit called, this is going to be the last split, keep split into 2 parts */
2084  if (nentries > INDEX_TUPLES_PER_PAGE && nentries <= 2 * INDEX_TUPLES_PER_PAGE)
2085  m = Max(m, nentries - INDEX_TUPLES_PER_PAGE);
2086 
2087  Assert(m >= 1);
2088 
2089  for (i = 0; i < (OffsetNumber)commonEntriesCount; i++)
2090  {
2091  BOX2DF *box = (BOX2DF *)DatumGetPointer(entryvec->vector[commonEntries[i].index].key);
2092 
2093  if (v->spl_nright > 0 && v->spl_nleft > 0)
2094  {
2095  /* Calculate delta between penalties of join "common entries" to different groups. */
2096  commonEntries[i].delta =
2097  fabsf(box2df_penalty(leftBox, box) - box2df_penalty(rightBox, box));
2098  }
2099  else
2100  {
2101  /* One group became degenerate during split - distribute boxes by their size */
2102  commonEntries[i].delta = -box2df_penalty_single(box);
2103  }
2104  }
2105 
2106  /* Sort "common entries" to distribute the most ambiguous entries first. */
2107  qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp);
2108 
2109  /* Distribute "common entries" between groups. */
2110  for (i = 0; i < (OffsetNumber)commonEntriesCount; i++)
2111  {
2112  bool place_right = true;
2113  BOX2DF *box = (BOX2DF *)DatumGetPointer(entryvec->vector[commonEntries[i].index].key);
2114 
2115  /* Recheck penalties. After we put undecided tuples to some side they're changed */
2116  float left_penalty = box2df_penalty(leftBox, box);
2117  float right_penalty = box2df_penalty(rightBox, box);
2118 
2119  /* Guarantee at least one tuple on each side to make Postgres happy */
2120  if (v->spl_nright == 0 && v->spl_nleft > 0)
2121  place_right = true;
2122  else if (v->spl_nleft == 0 && v->spl_nright > 0)
2123  place_right = false;
2124  /* Check if we have to place this entry in either group to achieve LIMIT_RATIO */
2125  else if (v->spl_nleft + (commonEntriesCount - i) <= m)
2126  place_right = false;
2127  else if (v->spl_nright + (commonEntriesCount - i) <= m)
2128  place_right = true;
2129  /* Otherwise select the group by smaller penalty */
2130  else if (left_penalty < right_penalty)
2131  place_right = false;
2132  else if (right_penalty < left_penalty)
2133  place_right = true;
2134  /* or just put tuple into smaller group */
2135  else if (v->spl_nleft < v->spl_nright)
2136  place_right = false;
2137  else
2138  place_right = true;
2139 
2140  if (place_right)
2141  PLACE_RIGHT(box, commonEntries[i].index);
2142  else
2143  PLACE_LEFT(box, commonEntries[i].index);
2144  }
2145  }
2146  v->spl_ldatum = PointerGetDatum(leftBox);
2147  v->spl_rdatum = PointerGetDatum(rightBox);
2148 
2149  POSTGIS_DEBUG(4, "[GIST] 'picksplit' completed");
2150 
2151  PG_RETURN_POINTER(v);
2152 }
2153 
2154 /*
2155 ** The BOX2DF key must be defined as a PostgreSQL type, even though it is only
2156 ** ever used internally. These no-op stubs are used to bind the type.
2157 */
2159 Datum box2df_in(PG_FUNCTION_ARGS)
2160 {
2161  ereport(ERROR,(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2162  errmsg("function box2df_in not implemented")));
2163  PG_RETURN_POINTER(NULL);
2164 }
2165 
2167 Datum box2df_out(PG_FUNCTION_ARGS)
2168 {
2169  BOX2DF *box = (BOX2DF *) PG_GETARG_POINTER(0);
2170  char *result = box2df_to_string(box);
2171  PG_RETURN_CSTRING(result);
2172 }
static uint8_t precision
Definition: cu_in_twkb.c:25
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
int gserialized_has_bbox(const GSERIALIZED *g)
Check if a GSERIALIZED has a bounding box without deserializing first.
Definition: gserialized.c:192
int32_t gserialized_get_srid(const GSERIALIZED *g)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
Definition: gserialized.c:155
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
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: gserialized.c:181
const float * gserialized_get_float_box_p(const GSERIALIZED *g, size_t *ndims)
Access to the float bounding box, if there is one.
Definition: gserialized.c:277
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
uint32_t gserialized_max_header_size(void)
Returns the size in bytes to read from toast to get the basic information from a geometry: GSERIALIZE...
Definition: gserialized.c:130
lwflags_t gserialized_get_lwflags(const GSERIALIZED *g)
Read the flags from a GSERIALIZED and return a standard lwflag integer.
Definition: gserialized.c:47
bool box2df_left(const BOX2DF *a, const BOX2DF *b)
static char * box2df_to_string(const BOX2DF *a)
void box2df_merge(BOX2DF *b_union, BOX2DF *b_new)
static void fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_sortsupport_2d(PG_FUNCTION_ARGS)
void box2df_set_empty(BOX2DF *a)
Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS)
Datum gserialized_contains_box2df_geom_2d(PG_FUNCTION_ARGS)
bool box2df_is_empty(const BOX2DF *a)
static int gserialized_datum_predicate_2d(Datum gs1, Datum gs2, box2df_predicate predicate)
Support function.
Datum gserialized_within_2d(PG_FUNCTION_ARGS)
PG_FUNCTION_INFO_V1(gserialized_contains_box2df_geom_2d)
Datum box2df_out(PG_FUNCTION_ARGS)
static Datum gserialized_gist_abbrev_convert_2d(Datum original, SortSupport ssup)
Datum gserialized_within_box2df_geom_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS)
Datum gserialized_left_2d(PG_FUNCTION_ARGS)
Datum gserialized_contains_2d(PG_FUNCTION_ARGS)
#define INDEX_TUPLES_PER_PAGE
Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS)
static int interval_cmp_upper(const void *i1, const void *i2)
Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS)
bool box2df_equals(const BOX2DF *a, const BOX2DF *b)
static uint64_t box2df_get_sortable_hash(const BOX2DF *b)
static int gserialized_gist_cmp_abbrev_2d(Datum x, Datum y, SortSupport ssup)
#define LIMIT_RATIO
bool box2df_contains(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS)
#define PLACE_RIGHT(box, off)
static int gserialized_gist_cmp_full_2d(Datum a, Datum b, SortSupport ssup)
Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS)
Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS)
Datum gserialized_contains_box2df_box2df_2d(PG_FUNCTION_ARGS)
Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS)
Datum gserialized_overlaps_box2df_geom_2d(PG_FUNCTION_ARGS)
Datum gserialized_within_box2df_box2df_2d(PG_FUNCTION_ARGS)
bool box2df_right(const BOX2DF *a, const BOX2DF *b)
bool box2df_overlaps(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_overleft_2d(PG_FUNCTION_ARGS)
void box2df_set_finite(BOX2DF *a)
static bool gserialized_gist_consistent_internal_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
Datum gserialized_overlaps_box2df_box2df_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
Datum gserialized_overright_2d(PG_FUNCTION_ARGS)
static bool box2df_within(const BOX2DF *a, const BOX2DF *b)
bool box2df_above(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_right_2d(PG_FUNCTION_ARGS)
static bool gserialized_gist_abbrev_abort_2d(int memtupcount, SortSupport ssup)
bool box2df_overbelow(const BOX2DF *a, const BOX2DF *b)
bool(* box2df_predicate)(const BOX2DF *a, const BOX2DF *b)
static double box2df_distance(const BOX2DF *a, const BOX2DF *b)
Calculate the box->box distance.
static double box2df_distance_leaf_centroid(const BOX2DF *a, const BOX2DF *b)
Calculate the centroid->centroid distance between the boxes.
Datum gserialized_below_2d(PG_FUNCTION_ARGS)
static float box2df_penalty_single(const BOX2DF *box)
Datum gserialized_overabove_2d(PG_FUNCTION_ARGS)
static void g_box_consider_split(ConsiderSplitContext *context, int dimNum, float rightLower, int minLeftCount, float leftUpper, int maxLeftCount)
int box2df_to_gbox_p(BOX2DF *a, GBOX *box)
static int interval_cmp_lower(const void *i1, const void *i2)
bool box2df_overright(const BOX2DF *a, const BOX2DF *b)
static double pt_distance(double ax, double ay, double bx, double by)
static float pack_float(const float value, const uint8_t realm)
static int common_entry_cmp(const void *i1, const void *i2)
Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS)
int gserialized_datum_get_box2df_p(Datum gsdatum, BOX2DF *box2df)
int gserialized_datum_get_internals_p(Datum gsdatum, GBOX *gbox, lwflags_t *flags, uint8_t *type, int32_t *srid)
Peak into a GSERIALIZED datum to find its bounding box and some other metadata.
int gserialized_datum_get_gbox_p(Datum gsdatum, GBOX *gbox)
Given a GSERIALIZED datum, as quickly as possible (peaking into the top of the memory) return the gbo...
static int box2df_from_gbox_p(GBOX *box, BOX2DF *a)
bool box2df_overabove(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS)
#define PLACE_LEFT(box, off)
void box2df_validate(BOX2DF *b)
static float box2df_penalty(const BOX2DF *b1, const BOX2DF *b2)
bool box2df_below(const BOX2DF *a, const BOX2DF *b)
static bool gserialized_gist_consistent_leaf_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
bool box2df_overleft(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_above_2d(PG_FUNCTION_ARGS)
static int gserialized_datum_predicate_box2df_geom_2d(const BOX2DF *br1, Datum gs2, box2df_predicate predicate)
BOX2DF * box2df_copy(BOX2DF *b)
Datum gserialized_same_2d(PG_FUNCTION_ARGS)
Datum box2df_in(PG_FUNCTION_ARGS)
static void adjustBox(BOX2DF *b, BOX2DF *addon)
static float non_negative(float val)
#define LW_FALSE
Definition: liblwgeom.h:94
#define LW_FAILURE
Definition: liblwgeom.h:96
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
Definition: liblwgeom.h:324
#define LW_SUCCESS
Definition: liblwgeom.h:97
uint16_t lwflags_t
Definition: liblwgeom.h:299
float next_float_up(double d)
Definition: lwgeom_api.c:74
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:93
float next_float_down(double d)
Definition: lwgeom_api.c:53
This library is the generic geometry handling section of PostGIS.
#define OUT_MAX_BYTES_DOUBLE
int lwprint_double(double d, int maxdd, char *buf)
Definition: lwprint.c:463
#define NAN
Definition: lwgeodetic.h:37
#define POSTGIS_FREE_IF_COPY_P(ptrsrc, ptrori)
Definition: lwinline.h:340
static uint64_t uint32_hilbert(uint32_t px, uint32_t py)
Definition: lwinline.h:256
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
int value
Definition: genraster.py:62
type
Definition: ovdump.py:42
double ymax
Definition: liblwgeom.h:357
double xmax
Definition: liblwgeom.h:355
double ymin
Definition: liblwgeom.h:356
double xmin
Definition: liblwgeom.h:354
uint32_t size
Definition: liblwgeom.h:444