PostGIS  2.3.8dev-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  *
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 
46 #include "../postgis_config.h"
47 
48 #include "liblwgeom.h" /* For standard geometry types. */
49 #include "lwgeom_pg.h" /* For debugging macros. */
50 #include "gserialized_gist.h" /* For utility functions. */
51 
52 /* Fall back to older finite() if necessary */
53 #ifndef HAVE_ISFINITE
54 # ifdef HAVE_GNU_ISFINITE
55 # define _GNU_SOURCE
56 # else
57 # define isfinite finite
58 # endif
59 #endif
60 
61 #include <float.h> /* For FLT_MAX */
62 
63 /*
64 ** When is a node split not so good? If more than 90% of the entries
65 ** end up in one of the children.
66 */
67 #define LIMIT_RATIO 0.1
68 
69 /*
70 ** 0 == don't use it
71 ** 1 == use it
72 */
73 #define KOROTKOV_SPLIT 1
74 
75 /*
76 ** For debugging
77 */
78 #if POSTGIS_DEBUG_LEVEL > 0
79 static int g2d_counter_leaf = 0;
80 static int g2d_counter_internal = 0;
81 #endif
82 
83 /*
84 ** GiST 2D key stubs
85 */
86 Datum box2df_out(PG_FUNCTION_ARGS);
87 Datum box2df_in(PG_FUNCTION_ARGS);
88 
89 /*
90 ** GiST 2D index function prototypes
91 */
92 Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS);
93 Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS);
94 Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS);
95 Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS);
96 Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS);
97 Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS);
98 Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS);
99 Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS);
100 
101 /*
102 ** GiST 2D operator prototypes
103 */
104 Datum gserialized_same_2d(PG_FUNCTION_ARGS);
105 Datum gserialized_within_2d(PG_FUNCTION_ARGS);
106 Datum gserialized_contains_2d(PG_FUNCTION_ARGS);
107 Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS);
108 Datum gserialized_left_2d(PG_FUNCTION_ARGS);
109 Datum gserialized_right_2d(PG_FUNCTION_ARGS);
110 Datum gserialized_above_2d(PG_FUNCTION_ARGS);
111 Datum gserialized_below_2d(PG_FUNCTION_ARGS);
112 Datum gserialized_overleft_2d(PG_FUNCTION_ARGS);
113 Datum gserialized_overright_2d(PG_FUNCTION_ARGS);
114 Datum gserialized_overabove_2d(PG_FUNCTION_ARGS);
115 Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS);
116 Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS);
117 Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS);
118 
119 #if POSTGIS_PGSQL_VERSION > 94
120 Datum gserialized_contains_box2df_geom_2d(PG_FUNCTION_ARGS);
121 Datum gserialized_contains_box2df_box2df_2d(PG_FUNCTION_ARGS);
122 Datum gserialized_within_box2df_geom_2d(PG_FUNCTION_ARGS);
123 Datum gserialized_within_box2df_box2df_2d(PG_FUNCTION_ARGS);
124 Datum gserialized_overlaps_box2df_geom_2d(PG_FUNCTION_ARGS);
125 Datum gserialized_overlaps_box2df_box2df_2d(PG_FUNCTION_ARGS);
126 #endif
127 
128 /*
129 ** true/false test function type
130 */
131 typedef bool (*box2df_predicate)(const BOX2DF *a, const BOX2DF *b);
132 
133 
134 static char* box2df_to_string(const BOX2DF *a)
135 {
136  char *rv = NULL;
137 
138  if ( a == NULL )
139  return pstrdup("<NULLPTR>");
140 
141  rv = palloc(128);
142  sprintf(rv, "BOX2DF(%.12g %.12g, %.12g %.12g)", a->xmin, a->ymin, a->xmax, a->ymax);
143  return rv;
144 }
145 
146 /* Allocate a new copy of BOX2DF */
147 static BOX2DF* box2df_copy(BOX2DF *b)
148 {
149  BOX2DF *c = (BOX2DF*)palloc(sizeof(BOX2DF));
150  memcpy((void*)c, (void*)b, sizeof(BOX2DF));
151  POSTGIS_DEBUGF(5, "copied box2df (%p) to box2df (%p)", b, c);
152  return c;
153 }
154 
155 
156 
157 /* Enlarge b_union to contain b_new. If b_new contains more
158  dimensions than b_union, expand b_union to contain those dimensions. */
159 static void box2df_merge(BOX2DF *b_union, BOX2DF *b_new)
160 {
161 
162  POSTGIS_DEBUGF(5, "merging %s with %s", box2df_to_string(b_union), box2df_to_string(b_new));
163  /* Adjust minimums */
164  if (b_union->xmin > b_new->xmin || isnan(b_union->xmin))
165  b_union->xmin = b_new->xmin;
166  if (b_union->ymin > b_new->ymin || isnan(b_union->ymin))
167  b_union->ymin = b_new->ymin;
168  /* Adjust maximums */
169  if (b_union->xmax < b_new->xmax || isnan(b_union->xmax))
170  b_union->xmax = b_new->xmax;
171  if (b_union->ymax < b_new->ymax || isnan(b_union->ymax))
172  b_union->ymax = b_new->ymax;
173 
174  POSTGIS_DEBUGF(5, "merge complete %s", box2df_to_string(b_union));
175  return;
176 }
177 
178 #if KOROTKOV_SPLIT < 1
179 static bool box2df_intersection(const BOX2DF *a, const BOX2DF *b, BOX2DF *n)
180 {
181  POSTGIS_DEBUGF(5, "calculating intersection of %s with %s", box2df_to_string(a), box2df_to_string(b));
182 
183  if( a == NULL || b == NULL || n == NULL )
184  return FALSE;
185 
186  n->xmax = Min(a->xmax, b->xmax);
187  n->ymax = Min(a->ymax, b->ymax);
188  n->xmin = Max(a->xmin, b->xmin);
189  n->ymin = Max(a->ymin, b->ymin);
190 
191  POSTGIS_DEBUGF(5, "intersection is %s", box2df_to_string(n));
192 
193  if ( (n->xmax < n->xmin) || (n->ymax < n->ymin) )
194  return FALSE;
195 
196  return TRUE;
197 }
198 #endif
199 
200 static float box2df_size(const BOX2DF *a)
201 {
202  float result;
203 
204  if ( a == NULL )
205  return (float)0.0;
206 
207  if ( (a->xmax <= a->xmin) || (a->ymax <= a->ymin) )
208  {
209  result = (float) 0.0;
210  }
211  else
212  {
213  result = (((double) a->xmax)-((double) a->xmin)) * (((double) a->ymax)-((double) a->ymin));
214  }
215 
216  return result;
217 }
218 
219 static float box2df_union_size(const BOX2DF *a, const BOX2DF *b)
220 {
221  float result;
222 
223  POSTGIS_DEBUG(5,"entered function");
224 
225  if ( a == NULL && b == NULL )
226  {
227  elog(ERROR, "box2df_union_size received two null arguments");
228  return 0.0;
229  }
230 
231  if ( a == NULL )
232  return box2df_size(b);
233 
234  if ( b == NULL )
235  return box2df_size(a);
236 
237  result = ((double)Max(a->xmax,b->xmax) - (double)Min(a->xmin,b->xmin)) *
238  ((double)Max(a->ymax,b->ymax) - (double)Min(a->ymin,b->ymin));
239 
240  POSTGIS_DEBUGF(5, "union size of %s and %s is %.8g", box2df_to_string(a), box2df_to_string(b), result);
241 
242  return result;
243 }
244 
245 
246 /* Convert a double-based GBOX into a float-based BOX2DF,
247  ensuring the float box is larger than the double box */
248 static inline int box2df_from_gbox_p(GBOX *box, BOX2DF *a)
249 {
250  a->xmin = next_float_down(box->xmin);
251  a->xmax = next_float_up(box->xmax);
252  a->ymin = next_float_down(box->ymin);
253  a->ymax = next_float_up(box->ymax);
254  return LW_SUCCESS;
255 }
256 
257 /***********************************************************************
258 ** BOX3DF tests for 2D index operators.
259 */
260 
261 /* Ensure all minimums are below maximums. */
262 static inline void box2df_validate(BOX2DF *b)
263 {
264  float tmp;
265  POSTGIS_DEBUGF(5,"validating box2df (%s)", box2df_to_string(b));
266  if ( b->xmax < b->xmin )
267  {
268  tmp = b->xmin;
269  b->xmin = b->xmax;
270  b->xmax = tmp;
271  }
272  if ( b->ymax < b->ymin )
273  {
274  tmp = b->ymin;
275  b->ymin = b->ymax;
276  b->ymax = tmp;
277  }
278  return;
279 }
280 
281 static bool box2df_overlaps(const BOX2DF *a, const BOX2DF *b)
282 {
283  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
284 
285  if ( (a->xmin > b->xmax) || (b->xmin > a->xmax) ||
286  (a->ymin > b->ymax) || (b->ymin > a->ymax) )
287  {
288  return FALSE;
289  }
290 
291  return TRUE;
292 }
293 
294 bool box2df_contains(const BOX2DF *a, const BOX2DF *b)
295 {
296  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
297 
298  if ( (a->xmin > b->xmin) || (a->xmax < b->xmax) ||
299  (a->ymin > b->ymin) || (a->ymax < b->ymax) )
300  {
301  return FALSE;
302  }
303 
304  return TRUE;
305 }
306 
307 static bool box2df_within(const BOX2DF *a, const BOX2DF *b)
308 {
309  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
310 
311  POSTGIS_DEBUG(5, "entered function");
312  return box2df_contains(b,a);
313 }
314 
315 static bool box2df_equals(const BOX2DF *a, const BOX2DF *b)
316 {
317  if ( a && b ) {
318  if ( (a->xmin != b->xmin) || (a->xmax != b->xmax) ||
319  (a->ymin != b->ymin) || (a->ymax != b->ymax) )
320  {
321  return FALSE;
322  }
323  return TRUE;
324  } else if ( a || b ) {
325  /* one empty, one not */
326  return FALSE;
327  } else {
328  /* both empty */
329  return TRUE;
330  }
331 }
332 
333 static bool box2df_overleft(const BOX2DF *a, const BOX2DF *b)
334 {
335  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
336 
337  /* a.xmax <= b.xmax */
338  return a->xmax <= b->xmax;
339 }
340 
341 static bool box2df_left(const BOX2DF *a, const BOX2DF *b)
342 {
343  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
344 
345  /* a.xmax < b.xmin */
346  return a->xmax < b->xmin;
347 }
348 
349 static bool box2df_right(const BOX2DF *a, const BOX2DF *b)
350 {
351  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
352 
353  /* a.xmin > b.xmax */
354  return a->xmin > b->xmax;
355 }
356 
357 static bool box2df_overright(const BOX2DF *a, const BOX2DF *b)
358 {
359  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
360 
361  /* a.xmin >= b.xmin */
362  return a->xmin >= b->xmin;
363 }
364 
365 static bool box2df_overbelow(const BOX2DF *a, const BOX2DF *b)
366 {
367  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
368 
369  /* a.ymax <= b.ymax */
370  return a->ymax <= b->ymax;
371 }
372 
373 static bool box2df_below(const BOX2DF *a, const BOX2DF *b)
374 {
375  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
376 
377  /* a.ymax < b.ymin */
378  return a->ymax < b->ymin;
379 }
380 
381 static bool box2df_above(const BOX2DF *a, const BOX2DF *b)
382 {
383  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
384 
385  /* a.ymin > b.ymax */
386  return a->ymin > b->ymax;
387 }
388 
389 static bool box2df_overabove(const BOX2DF *a, const BOX2DF *b)
390 {
391  if ( ! a || ! b ) return FALSE; /* TODO: might be smarter for EMPTY */
392 
393  /* a.ymin >= b.ymin */
394  return a->ymin >= b->ymin;
395 }
396 
400 static double box2df_distance_leaf_centroid(const BOX2DF *a, const BOX2DF *b)
401 {
402  /* The centroid->centroid distance between the boxes */
403  double a_x = (a->xmax + a->xmin) / 2.0;
404  double a_y = (a->ymax + a->ymin) / 2.0;
405  double b_x = (b->xmax + b->xmin) / 2.0;
406  double b_y = (b->ymax + b->ymin) / 2.0;
407 
408  /* This "distance" is only used for comparisons, */
409  return sqrt((a_x - b_x) * (a_x - b_x) + (a_y - b_y) * (a_y - b_y));
410 }
411 
412 #if POSTGIS_PGSQL_VERSION < 95
413 
417 static double box2df_distance_node_centroid(const BOX2DF *node, const BOX2DF *query)
418 {
419  BOX2DF q;
420  double qx, qy;
421  double d = 0.0;
422 
423  /* Turn query into point */
424  q.xmin = q.xmax = (query->xmin + query->xmax) / 2.0;
425  q.ymin = q.ymax = (query->ymin + query->ymax) / 2.0;
426  qx = q.xmin;
427  qy = q.ymin;
428 
429  /* Check for overlap */
430  if ( box2df_overlaps(node, &q) == LW_TRUE )
431  return 0.0;
432 
433  /* Above or below */
434  if ( qx >= node->xmin && qx <= node->xmax )
435  {
436  if( qy > node->ymax )
437  d = qy - node->ymax;
438  else if ( qy < node->ymin )
439  d = node->ymin - qy;
440  return d;
441  }
442  /* Left or right */
443  else if ( qy >= node->ymin && qy <= node->ymax )
444  {
445  if ( qx > node->xmax )
446  d = qx - node->xmax;
447  else if ( qx < node->xmin )
448  d = node->xmin - qx;
449  return d;
450  }
451  /* Corner quadrants */
452  else
453  {
454  /* below/left of xmin/ymin */
455  if ( qx < node->xmin && qy < node->ymin )
456  {
457  d = (node->xmin - qx) * (node->xmin - qx) +
458  (node->ymin - qy) * (node->ymin - qy);
459  }
460  /* above/left of xmin/ymax */
461  else if ( qx < node->xmin && qy > node->ymax )
462  {
463  d = (node->xmin - qx) * (node->xmin - qx) +
464  (node->ymax - qy) * (node->ymax - qy);
465  }
466  /* above/right of xmax/ymax */
467  else if ( qx > node->xmax && qy > node->ymax )
468  {
469  d = (node->xmax - qx) * (node->xmax - qx) +
470  (node->ymax - qy) * (node->ymax - qy);
471  }
472  /* below/right of xmax/ymin */
473  else if ( qx > node->xmin && qy < node->ymin )
474  {
475  d = (node->xmax - qx) * (node->xmax - qx) +
476  (node->ymin - qy) * (node->ymin - qy);
477  }
478  else
479  {
480  /*ERROR*/
481  elog(ERROR, "%s: reached unreachable code", __func__);
482  }
483  }
484 
485  return sqrt(d);
486 }
487 #endif
488 
489 /* Quick distance function */
490 static inline double pt_distance(double ax, double ay, double bx, double by)
491 {
492  return sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
493 }
494 
498 static double box2df_distance(const BOX2DF *a, const BOX2DF *b)
499 {
500  /* Check for overlap */
501  if ( box2df_overlaps(a, b) )
502  return 0.0;
503 
504  if ( box2df_left(a, b) )
505  {
506  if ( box2df_above(a, b) )
507  return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
508  if ( box2df_below(a, b) )
509  return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
510  else
511  return (double)b->xmin - (double)a->xmax;
512  }
513  if ( box2df_right(a, b) )
514  {
515  if ( box2df_above(a, b) )
516  return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
517  if ( box2df_below(a, b) )
518  return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
519  else
520  return (double)a->xmin - (double)b->xmax;
521  }
522  if ( box2df_above(a, b) )
523  {
524  if ( box2df_left(a, b) )
525  return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
526  if ( box2df_right(a, b) )
527  return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
528  else
529  return (double)a->ymin - (double)b->ymax;
530  }
531  if ( box2df_below(a, b) )
532  {
533  if ( box2df_left(a, b) )
534  return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
535  if ( box2df_right(a, b) )
536  return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
537  else
538  return (double)b->ymin - (double)a->ymax;
539  }
540 
541  return FLT_MAX;
542 }
543 
544 
551 int
552 gserialized_datum_get_box2df_p(Datum gsdatum, BOX2DF *box2df)
553 {
554  GSERIALIZED *gpart;
555  uint8_t flags;
556  int result = LW_SUCCESS;
557 
558  POSTGIS_DEBUG(4, "entered function");
559 
560  /*
561  ** The most info we need is the 8 bytes of serialized header plus the
562  ** of floats necessary to hold the bounding box.
563  */
564  if (VARATT_IS_EXTENDED(gsdatum))
565  {
566  gpart = (GSERIALIZED*)PG_DETOAST_DATUM_SLICE(gsdatum, 0, 8 + sizeof(BOX2DF));
567  }
568  else
569  {
570  gpart = (GSERIALIZED*)PG_DETOAST_DATUM(gsdatum);
571  }
572 
573  flags = gpart->flags;
574 
575  POSTGIS_DEBUGF(4, "got flags %d", gpart->flags);
576 
577  /* Do we even have a serialized bounding box? */
578  if ( FLAGS_GET_BBOX(flags) )
579  {
580  /* Yes! Copy it out into the box! */
581  POSTGIS_DEBUG(4, "copying box out of serialization");
582  memcpy(box2df, gpart->data, sizeof(BOX2DF));
583  result = LW_SUCCESS;
584  }
585  else
586  {
587  /* No, we need to calculate it from the full object. */
588  GBOX gbox;
589  GSERIALIZED *g = (GSERIALIZED*)PG_DETOAST_DATUM(gsdatum);
590 
591  gbox_init(&gbox);
592 
593  if (gserialized_get_gbox_p(g, &gbox) == LW_FAILURE)
594  {
595  POSTGIS_DEBUG(4, "could not calculate bbox, returning failure");
596  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
597  POSTGIS_FREE_IF_COPY_P(g, gsdatum);
598  return LW_FAILURE;
599  }
600  POSTGIS_FREE_IF_COPY_P(g, gsdatum);
601  result = box2df_from_gbox_p(&gbox, box2df);
602  }
603 
604  POSTGIS_FREE_IF_COPY_P(gpart, gsdatum);
605 
606  if ( result == LW_SUCCESS )
607  {
608  POSTGIS_DEBUGF(4, "got box2df %s", box2df_to_string(box2df));
609  }
610 
611  return result;
612 }
613 
614 
619 static int
620 gserialized_datum_predicate_2d(Datum gs1, Datum gs2, box2df_predicate predicate)
621 {
622  BOX2DF b1, b2, *br1=NULL, *br2=NULL;
623  POSTGIS_DEBUG(3, "entered function");
624 
625  if (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) br1 = &b1;
626  if (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) br2 = &b2;
627 
628  if ( predicate(br1, br2) )
629  {
630  POSTGIS_DEBUGF(3, "got boxes %s and %s", br1 ? box2df_to_string(&b1) : "(null)", br2 ? box2df_to_string(&b2) : "(null)");
631  return LW_TRUE;
632  }
633  return LW_FALSE;
634 }
635 
636 #if POSTGIS_PGSQL_VERSION > 94
637 static int
638 gserialized_datum_predicate_box2df_geom_2d(const BOX2DF *br1, Datum gs2, box2df_predicate predicate)
639 {
640  BOX2DF b2, *br2=NULL;
641  POSTGIS_DEBUG(3, "entered function");
642 
643  if (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) br2 = &b2;
644 
645  if ( predicate(br1, br2) )
646  {
647  POSTGIS_DEBUGF(3, "got boxes %s", br2 ? box2df_to_string(&b2) : "(null)");
648  return LW_TRUE;
649  }
650  return LW_FALSE;
651 }
652 
653 /***********************************************************************
654 * BRIN 2-D Index Operator Functions
655 */
656 
657 PG_FUNCTION_INFO_V1(gserialized_contains_box2df_geom_2d);
658 Datum gserialized_contains_box2df_geom_2d(PG_FUNCTION_ARGS)
659 {
660  POSTGIS_DEBUG(3, "entered function");
661  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_contains) == LW_TRUE )
662  PG_RETURN_BOOL(TRUE);
663 
664  PG_RETURN_BOOL(FALSE);
665 }
666 
667 PG_FUNCTION_INFO_V1(gserialized_contains_box2df_box2df_2d);
668 Datum gserialized_contains_box2df_box2df_2d(PG_FUNCTION_ARGS)
669 {
670  if ( box2df_contains((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
671  PG_RETURN_BOOL(TRUE);
672 
673  PG_RETURN_BOOL(FALSE);
674 }
675 
676 PG_FUNCTION_INFO_V1(gserialized_within_box2df_geom_2d);
677 Datum gserialized_within_box2df_geom_2d(PG_FUNCTION_ARGS)
678 {
679  POSTGIS_DEBUG(3, "entered function");
680  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_within) == LW_TRUE )
681  PG_RETURN_BOOL(TRUE);
682 
683  PG_RETURN_BOOL(FALSE);
684 }
685 
686 PG_FUNCTION_INFO_V1(gserialized_within_box2df_box2df_2d);
687 Datum gserialized_within_box2df_box2df_2d(PG_FUNCTION_ARGS)
688 {
689  if ( box2df_within((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
690  PG_RETURN_BOOL(TRUE);
691 
692  PG_RETURN_BOOL(FALSE);
693 }
694 
695 PG_FUNCTION_INFO_V1(gserialized_overlaps_box2df_geom_2d);
696 Datum gserialized_overlaps_box2df_geom_2d(PG_FUNCTION_ARGS)
697 {
698  POSTGIS_DEBUG(3, "entered function");
699  if ( gserialized_datum_predicate_box2df_geom_2d((BOX2DF*)PG_GETARG_POINTER(0), PG_GETARG_DATUM(1), box2df_overlaps) == LW_TRUE )
700  PG_RETURN_BOOL(TRUE);
701 
702  PG_RETURN_BOOL(FALSE);
703 }
704 
705 PG_FUNCTION_INFO_V1(gserialized_overlaps_box2df_box2df_2d);
706 Datum gserialized_overlaps_box2df_box2df_2d(PG_FUNCTION_ARGS)
707 {
708  if ( box2df_overlaps((BOX2DF *)PG_GETARG_POINTER(0), (BOX2DF *)PG_GETARG_POINTER(1)))
709  PG_RETURN_BOOL(TRUE);
710 
711  PG_RETURN_BOOL(FALSE);
712 }
713 #endif
714 
715 /***********************************************************************
716 * GiST 2-D Index Operator Functions
717 */
718 
720 Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS)
721 {
722  BOX2DF b1, b2;
723  Datum gs1 = PG_GETARG_DATUM(0);
724  Datum gs2 = PG_GETARG_DATUM(1);
725 
726  POSTGIS_DEBUG(3, "entered function");
727 
728  /* Must be able to build box for each argument (ie, not empty geometry). */
729  if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
731  {
732  double distance = box2df_distance_leaf_centroid(&b1, &b2);
733  POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
734  PG_RETURN_FLOAT8(distance);
735  }
736  PG_RETURN_FLOAT8(FLT_MAX);
737 }
738 
740 Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS)
741 {
742  BOX2DF b1, b2;
743  Datum gs1 = PG_GETARG_DATUM(0);
744  Datum gs2 = PG_GETARG_DATUM(1);
745 
746  POSTGIS_DEBUG(3, "entered function");
747 
748  /* Must be able to build box for each argument (ie, not empty geometry). */
749  if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
751  {
752  double distance = box2df_distance(&b1, &b2);
753  POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
754  PG_RETURN_FLOAT8(distance);
755  }
756  PG_RETURN_FLOAT8(FLT_MAX);
757 }
758 
760 Datum gserialized_same_2d(PG_FUNCTION_ARGS)
761 {
762  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_equals) == LW_TRUE )
763  PG_RETURN_BOOL(TRUE);
764 
765  PG_RETURN_BOOL(FALSE);
766 }
767 
769 Datum gserialized_within_2d(PG_FUNCTION_ARGS)
770 {
771  POSTGIS_DEBUG(3, "entered function");
772  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_within) == LW_TRUE )
773  PG_RETURN_BOOL(TRUE);
774 
775  PG_RETURN_BOOL(FALSE);
776 }
777 
779 Datum gserialized_contains_2d(PG_FUNCTION_ARGS)
780 {
781  POSTGIS_DEBUG(3, "entered function");
782  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_contains) == LW_TRUE )
783  PG_RETURN_BOOL(TRUE);
784 
785  PG_RETURN_BOOL(FALSE);
786 }
787 
789 Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS)
790 {
791  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overlaps) == LW_TRUE )
792  PG_RETURN_BOOL(TRUE);
793 
794  PG_RETURN_BOOL(FALSE);
795 }
796 
798 Datum gserialized_left_2d(PG_FUNCTION_ARGS)
799 {
800  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_left) == LW_TRUE )
801  PG_RETURN_BOOL(TRUE);
802 
803  PG_RETURN_BOOL(FALSE);
804 }
805 
807 Datum gserialized_right_2d(PG_FUNCTION_ARGS)
808 {
809  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_right) == LW_TRUE )
810  PG_RETURN_BOOL(TRUE);
811 
812  PG_RETURN_BOOL(FALSE);
813 }
814 
816 Datum gserialized_above_2d(PG_FUNCTION_ARGS)
817 {
818  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_above) == LW_TRUE )
819  PG_RETURN_BOOL(TRUE);
820 
821  PG_RETURN_BOOL(FALSE);
822 }
823 
825 Datum gserialized_below_2d(PG_FUNCTION_ARGS)
826 {
827  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_below) == LW_TRUE )
828  PG_RETURN_BOOL(TRUE);
829 
830  PG_RETURN_BOOL(FALSE);
831 }
832 
834 Datum gserialized_overleft_2d(PG_FUNCTION_ARGS)
835 {
836  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overleft) == LW_TRUE )
837  PG_RETURN_BOOL(TRUE);
838 
839  PG_RETURN_BOOL(FALSE);
840 }
841 
843 Datum gserialized_overright_2d(PG_FUNCTION_ARGS)
844 {
845  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overright) == LW_TRUE )
846  PG_RETURN_BOOL(TRUE);
847 
848  PG_RETURN_BOOL(FALSE);
849 }
850 
852 Datum gserialized_overabove_2d(PG_FUNCTION_ARGS)
853 {
854  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overabove) == LW_TRUE )
855  PG_RETURN_BOOL(TRUE);
856 
857  PG_RETURN_BOOL(FALSE);
858 }
859 
861 Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS)
862 {
863  if ( gserialized_datum_predicate_2d(PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), box2df_overbelow) == LW_TRUE )
864  PG_RETURN_BOOL(TRUE);
865 
866  PG_RETURN_BOOL(FALSE);
867 }
868 
869 
870 /***********************************************************************
871 * GiST Index Support Functions
872 */
873 
874 /*
875 ** GiST support function. Given a geography, return a "compressed"
876 ** version. In this case, we convert the geography into a geocentric
877 ** bounding box. If the geography already has the box embedded in it
878 ** we pull that out and hand it back.
879 */
881 Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS)
882 {
883  GISTENTRY *entry_in = (GISTENTRY*)PG_GETARG_POINTER(0);
884  GISTENTRY *entry_out = NULL;
885  BOX2DF bbox_out;
886  int result = LW_SUCCESS;
887 
888  POSTGIS_DEBUG(4, "[GIST] 'compress' function called");
889 
890  /*
891  ** Not a leaf key? There's nothing to do.
892  ** Return the input unchanged.
893  */
894  if ( ! entry_in->leafkey )
895  {
896  POSTGIS_DEBUG(4, "[GIST] non-leafkey entry, returning input unaltered");
897  PG_RETURN_POINTER(entry_in);
898  }
899 
900  POSTGIS_DEBUG(4, "[GIST] processing leafkey input");
901  entry_out = palloc(sizeof(GISTENTRY));
902 
903  /*
904  ** Null key? Make a copy of the input entry and
905  ** return.
906  */
907  if ( DatumGetPointer(entry_in->key) == NULL )
908  {
909  POSTGIS_DEBUG(4, "[GIST] leafkey is null");
910  gistentryinit(*entry_out, (Datum) 0, entry_in->rel,
911  entry_in->page, entry_in->offset, FALSE);
912  POSTGIS_DEBUG(4, "[GIST] returning copy of input");
913  PG_RETURN_POINTER(entry_out);
914  }
915 
916  /* Extract our index key from the GiST entry. */
917  result = gserialized_datum_get_box2df_p(entry_in->key, &bbox_out);
918 
919  /* Is the bounding box valid (non-empty, non-infinite)? If not, return input uncompressed. */
920  if ( result == LW_FAILURE )
921  {
922  POSTGIS_DEBUG(4, "[GIST] empty geometry!");
923  PG_RETURN_POINTER(entry_in);
924  }
925 
926  POSTGIS_DEBUGF(4, "[GIST] got entry_in->key: %s", box2df_to_string(&bbox_out));
927 
928  /* Check all the dimensions for finite values */
929  if ( ! isfinite(bbox_out.xmax) || ! isfinite(bbox_out.xmin) ||
930  ! isfinite(bbox_out.ymax) || ! isfinite(bbox_out.ymin) )
931  {
932  POSTGIS_DEBUG(4, "[GIST] infinite geometry!");
933  PG_RETURN_POINTER(entry_in);
934  }
935 
936  /* Enure bounding box has minimums below maximums. */
937  box2df_validate(&bbox_out);
938 
939  /* Prepare GISTENTRY for return. */
940  gistentryinit(*entry_out, PointerGetDatum(box2df_copy(&bbox_out)),
941  entry_in->rel, entry_in->page, entry_in->offset, FALSE);
942 
943  /* Return GISTENTRY. */
944  POSTGIS_DEBUG(4, "[GIST] 'compress' function complete");
945  PG_RETURN_POINTER(entry_out);
946 }
947 
948 
949 /*
950 ** GiST support function.
951 ** Decompress an entry. Unused for geography, so we return.
952 */
954 Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS)
955 {
956  POSTGIS_DEBUG(5, "[GIST] 'decompress' function called");
957  /* We don't decompress. Just return the input. */
958  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
959 }
960 
961 
962 
963 /*
964 ** GiST support function. Called from gserialized_gist_consistent below.
965 */
966 static inline bool gserialized_gist_consistent_leaf_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
967 {
968  bool retval;
969 
970  POSTGIS_DEBUGF(4, "[GIST] leaf consistent, strategy [%d], count[%i]",
971  strategy, g2d_counter_leaf++);
972 
973  switch (strategy)
974  {
975 
976  /* Basic overlaps */
977  case RTOverlapStrategyNumber:
978  retval = (bool) box2df_overlaps(key, query);
979  break;
980  case RTSameStrategyNumber:
981  retval = (bool) box2df_equals(key, query);
982  break;
983  case RTContainsStrategyNumber:
984  case RTOldContainsStrategyNumber:
985  retval = (bool) box2df_contains(key, query);
986  break;
987  case RTContainedByStrategyNumber:
988  case RTOldContainedByStrategyNumber:
989  retval = (bool) box2df_contains(query, key);
990  break;
991 
992  /* To one side */
993  case RTAboveStrategyNumber:
994  retval = (bool) box2df_above(key, query);
995  break;
996  case RTBelowStrategyNumber:
997  retval = (bool) box2df_below(key, query);
998  break;
999  case RTRightStrategyNumber:
1000  retval = (bool) box2df_right(key, query);
1001  break;
1002  case RTLeftStrategyNumber:
1003  retval = (bool) box2df_left(key, query);
1004  break;
1005 
1006  /* Overlapping to one side */
1007  case RTOverAboveStrategyNumber:
1008  retval = (bool) box2df_overabove(key, query);
1009  break;
1010  case RTOverBelowStrategyNumber:
1011  retval = (bool) box2df_overbelow(key, query);
1012  break;
1013  case RTOverRightStrategyNumber:
1014  retval = (bool) box2df_overright(key, query);
1015  break;
1016  case RTOverLeftStrategyNumber:
1017  retval = (bool) box2df_overleft(key, query);
1018  break;
1019 
1020  default:
1021  retval = FALSE;
1022  }
1023 
1024  return (retval);
1025 }
1026 
1027 /*
1028 ** GiST support function. Called from gserialized_gist_consistent below.
1029 */
1030 static inline bool gserialized_gist_consistent_internal_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
1031 {
1032  bool retval;
1033 
1034  POSTGIS_DEBUGF(4, "[GIST] internal consistent, strategy [%d], count[%i], query[%s], key[%s]",
1035  strategy, g2d_counter_internal++, box2df_to_string(query), box2df_to_string(key) );
1036 
1037  switch (strategy)
1038  {
1039 
1040  /* Basic overlaps */
1041  case RTOverlapStrategyNumber:
1042  retval = (bool) box2df_overlaps(key, query);
1043  break;
1044  case RTSameStrategyNumber:
1045  case RTContainsStrategyNumber:
1046  case RTOldContainsStrategyNumber:
1047  retval = (bool) box2df_contains(key, query);
1048  break;
1049  case RTContainedByStrategyNumber:
1050  case RTOldContainedByStrategyNumber:
1051  retval = (bool) box2df_overlaps(key, query);
1052  break;
1053 
1054  /* To one side */
1055  case RTAboveStrategyNumber:
1056  retval = (bool)(!box2df_overbelow(key, query));
1057  break;
1058  case RTBelowStrategyNumber:
1059  retval = (bool)(!box2df_overabove(key, query));
1060  break;
1061  case RTRightStrategyNumber:
1062  retval = (bool)(!box2df_overleft(key, query));
1063  break;
1064  case RTLeftStrategyNumber:
1065  retval = (bool)(!box2df_overright(key, query));
1066  break;
1067 
1068  /* Overlapping to one side */
1069  case RTOverAboveStrategyNumber:
1070  retval = (bool)(!box2df_below(key, query));
1071  break;
1072  case RTOverBelowStrategyNumber:
1073  retval = (bool)(!box2df_above(key, query));
1074  break;
1075  case RTOverRightStrategyNumber:
1076  retval = (bool)(!box2df_left(key, query));
1077  break;
1078  case RTOverLeftStrategyNumber:
1079  retval = (bool)(!box2df_right(key, query));
1080  break;
1081 
1082  default:
1083  retval = FALSE;
1084  }
1085 
1086  return (retval);
1087 }
1088 
1089 /*
1090 ** GiST support function. Take in a query and an entry and see what the
1091 ** relationship is, based on the query strategy.
1092 */
1094 Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS)
1095 {
1096  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1097  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1098  bool result;
1099  BOX2DF query_gbox_index;
1100 
1101  /* PostgreSQL 8.4 and later require the RECHECK flag to be set here,
1102  rather than being supplied as part of the operator class definition */
1103  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1104 
1105  /* We set recheck to false to avoid repeatedly pulling every "possibly matched" geometry
1106  out during index scans. For cases when the geometries are large, rechecking
1107  can make things twice as slow. */
1108  *recheck = false;
1109 
1110  POSTGIS_DEBUG(4, "[GIST] 'consistent' function called");
1111 
1112  /* Quick sanity check on query argument. */
1113  if ( DatumGetPointer(PG_GETARG_DATUM(1)) == NULL )
1114  {
1115  POSTGIS_DEBUG(4, "[GIST] null query pointer (!?!), returning false");
1116  PG_RETURN_BOOL(FALSE); /* NULL query! This is screwy! */
1117  }
1118 
1119  /* Quick sanity check on entry key. */
1120  if ( DatumGetPointer(entry->key) == NULL )
1121  {
1122  POSTGIS_DEBUG(4, "[GIST] null index entry, returning false");
1123  PG_RETURN_BOOL(FALSE); /* NULL entry! */
1124  }
1125 
1126  /* Null box should never make this far. */
1127  if ( gserialized_datum_get_box2df_p(PG_GETARG_DATUM(1), &query_gbox_index) == LW_FAILURE )
1128  {
1129  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1130  PG_RETURN_BOOL(FALSE);
1131  }
1132 
1133  /* Treat leaf node tests different from internal nodes */
1134  if (GIST_LEAF(entry))
1135  {
1137  (BOX2DF*)DatumGetPointer(entry->key),
1138  &query_gbox_index, strategy);
1139  }
1140  else
1141  {
1143  (BOX2DF*)DatumGetPointer(entry->key),
1144  &query_gbox_index, strategy);
1145  }
1146 
1147  PG_RETURN_BOOL(result);
1148 }
1149 
1150 
1151 /*
1152 ** GiST support function. Take in a query and an entry and return the "distance"
1153 ** between them.
1154 **
1155 ** Given an index entry p and a query value q, this function determines the
1156 ** index entry's "distance" from the query value. This function must be
1157 ** supplied if the operator class contains any ordering operators. A query
1158 ** using the ordering operator will be implemented by returning index entries
1159 ** with the smallest "distance" values first, so the results must be consistent
1160 ** with the operator's semantics. For a leaf index entry the result just
1161 ** represents the distance to the index entry; for an internal tree node, the
1162 ** result must be the smallest distance that any child entry could have.
1163 **
1164 ** Strategy 13 = true distance tests <->
1165 ** Strategy 14 = box-based distance tests <#>
1166 */
1168 Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
1169 {
1170  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1171  BOX2DF query_box;
1172  BOX2DF *entry_box;
1173  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1174  double distance;
1175 #if POSTGIS_PGSQL_VERSION >= 95
1176  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1177 #endif
1178 
1179  POSTGIS_DEBUG(4, "[GIST] 'distance' function called");
1180 
1181  /* We are using '13' as the gist true-distance <-> strategy number
1182  * and '14' as the gist distance-between-boxes <#> strategy number */
1183  if ( strategy != 13 && strategy != 14 ) {
1184  elog(ERROR, "unrecognized strategy number: %d", strategy);
1185  PG_RETURN_FLOAT8(FLT_MAX);
1186  }
1187 
1188  /* Null box should never make this far. */
1189  if ( gserialized_datum_get_box2df_p(PG_GETARG_DATUM(1), &query_box) == LW_FAILURE )
1190  {
1191  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1192  PG_RETURN_FLOAT8(FLT_MAX);
1193  }
1194 
1195  /* Get the entry box */
1196  entry_box = (BOX2DF*)DatumGetPointer(entry->key);
1197 
1198 #if POSTGIS_PGSQL_VERSION >= 95
1199 
1200  /* Box-style distance test */
1201  if ( strategy == 14 ) /* operator <#> */
1202  {
1203  distance = box2df_distance(entry_box, &query_box);
1204  }
1205  /* True distance test (formerly centroid distance) */
1206  else if ( strategy == 13 ) /* operator <-> */
1207  {
1208  /* In all cases, since we only have keys (boxes) we'll return */
1209  /* the minimum possible distance, which is the box2df_distance */
1210  /* and let the recheck sort things out in the case of leaves */
1211  distance = box2df_distance(entry_box, &query_box);
1212 
1213  if (GIST_LEAF(entry))
1214  *recheck = true;
1215  }
1216  else
1217  {
1218  elog(ERROR, "%s: reached unreachable code", __func__);
1219  PG_RETURN_NULL();
1220  }
1221 #else
1222  /* Box-style distance test */
1223  if ( strategy == 14 )
1224  {
1225  distance = (double)box2df_distance(entry_box, &query_box);
1226  PG_RETURN_FLOAT8(distance);
1227  }
1228 
1229  /* Treat leaf node tests different from internal nodes */
1230  if (GIST_LEAF(entry))
1231  {
1232  /* Calculate distance to leaves */
1233  distance = (double)box2df_distance_leaf_centroid(entry_box, &query_box);
1234  }
1235  else
1236  {
1237  /* Calculate distance for internal nodes */
1238  distance = (double)box2df_distance_node_centroid(entry_box, &query_box);
1239  }
1240 #endif
1241 
1242  PG_RETURN_FLOAT8(distance);
1243 }
1244 
1245 /*
1246 ** GiST support function. Calculate the "penalty" cost of adding this entry into an existing entry.
1247 ** Calculate the change in volume of the old entry once the new entry is added.
1248 ** TODO: Re-evaluate this in light of R*Tree penalty approaches.
1249 */
1251 Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS)
1252 {
1253  GISTENTRY *origentry = (GISTENTRY*) PG_GETARG_POINTER(0);
1254  GISTENTRY *newentry = (GISTENTRY*) PG_GETARG_POINTER(1);
1255  float *result = (float*) PG_GETARG_POINTER(2);
1256  BOX2DF *gbox_index_orig, *gbox_index_new;
1257  float size_union, size_orig;
1258 
1259  POSTGIS_DEBUG(4, "[GIST] 'penalty' function called");
1260 
1261  gbox_index_orig = (BOX2DF*)DatumGetPointer(origentry->key);
1262  gbox_index_new = (BOX2DF*)DatumGetPointer(newentry->key);
1263 
1264  /* Drop out if we're dealing with null inputs. Shouldn't happen. */
1265  if ( (gbox_index_orig == NULL) && (gbox_index_new == NULL) )
1266  {
1267  POSTGIS_DEBUG(4, "[GIST] both inputs NULL! returning penalty of zero");
1268  *result = 0.0;
1269  PG_RETURN_FLOAT8(*result);
1270  }
1271 
1272  /* Calculate the size difference of the boxes. */
1273  size_union = box2df_union_size(gbox_index_orig, gbox_index_new);
1274  size_orig = box2df_size(gbox_index_orig);
1275  *result = size_union - size_orig;
1276 
1277  POSTGIS_DEBUGF(4, "[GIST] 'penalty', union size (%.12f), original size (%.12f), penalty (%.12f)", size_union, size_orig, *result);
1278 
1279  PG_RETURN_POINTER(result);
1280 }
1281 
1282 /*
1283 ** GiST support function. Merge all the boxes in a page into one master box.
1284 */
1286 Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS)
1287 {
1288  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
1289  int *sizep = (int *) PG_GETARG_POINTER(1); /* Size of the return value */
1290  int numranges, i;
1291  BOX2DF *box_cur, *box_union;
1292 
1293  POSTGIS_DEBUG(4, "[GIST] 'union' function called");
1294 
1295  numranges = entryvec->n;
1296 
1297  box_cur = (BOX2DF*) DatumGetPointer(entryvec->vector[0].key);
1298 
1299  box_union = box2df_copy(box_cur);
1300 
1301  for ( i = 1; i < numranges; i++ )
1302  {
1303  box_cur = (BOX2DF*) DatumGetPointer(entryvec->vector[i].key);
1304  box2df_merge(box_union, box_cur);
1305  }
1306 
1307  *sizep = sizeof(BOX2DF);
1308 
1309  POSTGIS_DEBUGF(4, "[GIST] 'union', numranges(%i), pageunion %s", numranges, box2df_to_string(box_union));
1310 
1311  PG_RETURN_POINTER(box_union);
1312 
1313 }
1314 
1315 /*
1316 ** GiST support function. Test equality of keys.
1317 */
1319 Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS)
1320 {
1321  BOX2DF *b1 = (BOX2DF*)PG_GETARG_POINTER(0);
1322  BOX2DF *b2 = (BOX2DF*)PG_GETARG_POINTER(1);
1323  bool *result = (bool*)PG_GETARG_POINTER(2);
1324 
1325  POSTGIS_DEBUG(4, "[GIST] 'same' function called");
1326 
1327  *result = box2df_equals(b1, b2);
1328 
1329  PG_RETURN_POINTER(result);
1330 }
1331 
1332 #if KOROTKOV_SPLIT > 0
1333 /*
1334  * Adjust BOX2DF b boundaries with insertion of addon.
1335  */
1336 static void
1337 adjustBox(BOX2DF *b, BOX2DF *addon)
1338 {
1339  if (b->xmax < addon->xmax || isnan(b->xmax))
1340  b->xmax = addon->xmax;
1341  if (b->xmin > addon->xmin || isnan(b->xmin))
1342  b->xmin = addon->xmin;
1343  if (b->ymax < addon->ymax || isnan(b->ymax))
1344  b->ymax = addon->ymax;
1345  if (b->ymin > addon->ymin || isnan(b->ymin))
1346  b->ymin = addon->ymin;
1347 }
1348 
1349 /*
1350  * Trivial split: half of entries will be placed on one page
1351  * and another half - to another
1352  */
1353 static void
1354 fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
1355 {
1356  OffsetNumber i,
1357  maxoff;
1358  BOX2DF *unionL = NULL,
1359  *unionR = NULL;
1360  int nbytes;
1361 
1362  maxoff = entryvec->n - 1;
1363 
1364  nbytes = (maxoff + 2) * sizeof(OffsetNumber);
1365  v->spl_left = (OffsetNumber *) palloc(nbytes);
1366  v->spl_right = (OffsetNumber *) palloc(nbytes);
1367  v->spl_nleft = v->spl_nright = 0;
1368 
1369  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1370  {
1371  BOX2DF *cur = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
1372 
1373  if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
1374  {
1375  v->spl_left[v->spl_nleft] = i;
1376  if (unionL == NULL)
1377  {
1378  unionL = (BOX2DF *) palloc(sizeof(BOX2DF));
1379  *unionL = *cur;
1380  }
1381  else
1382  adjustBox(unionL, cur);
1383 
1384  v->spl_nleft++;
1385  }
1386  else
1387  {
1388  v->spl_right[v->spl_nright] = i;
1389  if (unionR == NULL)
1390  {
1391  unionR = (BOX2DF *) palloc(sizeof(BOX2DF));
1392  *unionR = *cur;
1393  }
1394  else
1395  adjustBox(unionR, cur);
1396 
1397  v->spl_nright++;
1398  }
1399  }
1400 
1401  if (v->spl_ldatum_exists)
1402  adjustBox(unionL, (BOX2DF *) DatumGetPointer(v->spl_ldatum));
1403  v->spl_ldatum = BoxPGetDatum(unionL);
1404 
1405  if (v->spl_rdatum_exists)
1406  adjustBox(unionR, (BOX2DF *) DatumGetPointer(v->spl_rdatum));
1407  v->spl_rdatum = BoxPGetDatum(unionR);
1408 
1409  v->spl_ldatum_exists = v->spl_rdatum_exists = false;
1410 }
1411 
1412 /*
1413  * Represents information about an entry that can be placed to either group
1414  * without affecting overlap over selected axis ("common entry").
1415  */
1416 typedef struct
1417 {
1418  /* Index of entry in the initial array */
1419  int index;
1420  /* Delta between penalties of entry insertion into different groups */
1421  float delta;
1422 } CommonEntry;
1423 
1424 /*
1425  * Context for g_box_consider_split. Contains information about currently
1426  * selected split and some general information.
1427  */
1428 typedef struct
1429 {
1430  int entriesCount; /* total number of entries being split */
1431  BOX2DF boundingBox; /* minimum bounding box across all entries */
1432 
1433  /* Information about currently selected split follows */
1434 
1435  bool first; /* true if no split was selected yet */
1436 
1437  float leftUpper; /* upper bound of left interval */
1438  float rightLower; /* lower bound of right interval */
1439 
1440  float4 ratio;
1441  float4 overlap;
1442  int dim; /* axis of this split */
1443  float range; /* width of general MBR projection to the
1444  * selected axis */
1446 
1447 /*
1448  * Interval represents projection of box to axis.
1449  */
1450 typedef struct
1451 {
1452  float lower,
1453  upper;
1454 } SplitInterval;
1455 
1456 /*
1457  * Interval comparison function by lower bound of the interval;
1458  */
1459 static int
1460 interval_cmp_lower(const void *i1, const void *i2)
1461 {
1462  float lower1 = ((const SplitInterval *) i1)->lower,
1463  lower2 = ((const SplitInterval *) i2)->lower;
1464 
1465  if (isnan(lower1))
1466  {
1467  if (isnan(lower2))
1468  return 0;
1469  else
1470  return 1;
1471  }
1472  else if (isnan(lower2))
1473  {
1474  return -1;
1475  }
1476  else
1477  {
1478  if (lower1 < lower2)
1479  return -1;
1480  else if (lower1 > lower2)
1481  return 1;
1482  else
1483  return 0;
1484  }
1485 }
1486 
1487 
1488 /*
1489  * Interval comparison function by upper bound of the interval;
1490  */
1491 static int
1492 interval_cmp_upper(const void *i1, const void *i2)
1493 {
1494  float upper1 = ((const SplitInterval *) i1)->upper,
1495  upper2 = ((const SplitInterval *) i2)->upper;
1496 
1497  if (isnan(upper1))
1498  {
1499  if (isnan(upper2))
1500  return 0;
1501  else
1502  return -1;
1503  }
1504  else if (isnan(upper2))
1505  {
1506  return 1;
1507  }
1508  else
1509  {
1510  if (upper1 < upper2)
1511  return -1;
1512  else if (upper1 > upper2)
1513  return 1;
1514  else
1515  return 0;
1516  }
1517 }
1518 
1519 /*
1520  * Replace negative value with zero.
1521  */
1522 static inline float
1523 non_negative(float val)
1524 {
1525  if (val >= 0.0f)
1526  return val;
1527  else
1528  return 0.0f;
1529 }
1530 
1531 /*
1532  * Consider replacement of currently selected split with the better one.
1533  */
1534 static inline void
1536  float rightLower, int minLeftCount,
1537  float leftUpper, int maxLeftCount)
1538 {
1539  int leftCount,
1540  rightCount;
1541  float4 ratio,
1542  overlap;
1543  float range;
1544 
1545  POSTGIS_DEBUGF(5, "consider split: dimNum = %d, rightLower = %f, "
1546  "minLeftCount = %d, leftUpper = %f, maxLeftCount = %d ",
1547  dimNum, rightLower, minLeftCount, leftUpper, maxLeftCount);
1548 
1549  /*
1550  * Calculate entries distribution ratio assuming most uniform distribution
1551  * of common entries.
1552  */
1553  if (minLeftCount >= (context->entriesCount + 1) / 2)
1554  {
1555  leftCount = minLeftCount;
1556  }
1557  else
1558  {
1559  if (maxLeftCount <= context->entriesCount / 2)
1560  leftCount = maxLeftCount;
1561  else
1562  leftCount = context->entriesCount / 2;
1563  }
1564  rightCount = context->entriesCount - leftCount;
1565 
1566  /*
1567  * Ratio of split - quotient between size of lesser group and total
1568  * entries count.
1569  */
1570  ratio = ((float4) Min(leftCount, rightCount)) /
1571  ((float4) context->entriesCount);
1572 
1573  if (ratio > LIMIT_RATIO)
1574  {
1575  bool selectthis = false;
1576 
1577  /*
1578  * The ratio is acceptable, so compare current split with previously
1579  * selected one. Between splits of one dimension we search for minimal
1580  * overlap (allowing negative values) and minimal ration (between same
1581  * overlaps. We switch dimension if find less overlap (non-negative)
1582  * or less range with same overlap.
1583  */
1584  if (dimNum == 0)
1585  range = context->boundingBox.xmax - context->boundingBox.xmin;
1586  else
1587  range = context->boundingBox.ymax - context->boundingBox.ymin;
1588 
1589  overlap = (leftUpper - rightLower) / range;
1590 
1591  /* If there is no previous selection, select this */
1592  if (context->first)
1593  selectthis = true;
1594  else if (context->dim == dimNum)
1595  {
1596  /*
1597  * Within the same dimension, choose the new split if it has a
1598  * smaller overlap, or same overlap but better ratio.
1599  */
1600  if (overlap < context->overlap ||
1601  (overlap == context->overlap && ratio > context->ratio))
1602  selectthis = true;
1603  }
1604  else
1605  {
1606  /*
1607  * Across dimensions, choose the new split if it has a smaller
1608  * *non-negative* overlap, or same *non-negative* overlap but
1609  * bigger range. This condition differs from the one described in
1610  * the article. On the datasets where leaf MBRs don't overlap
1611  * themselves, non-overlapping splits (i.e. splits which have zero
1612  * *non-negative* overlap) are frequently possible. In this case
1613  * splits tends to be along one dimension, because most distant
1614  * non-overlapping splits (i.e. having lowest negative overlap)
1615  * appears to be in the same dimension as in the previous split.
1616  * Therefore MBRs appear to be very prolonged along another
1617  * dimension, which leads to bad search performance. Using range
1618  * as the second split criteria makes MBRs more quadratic. Using
1619  * *non-negative* overlap instead of overlap as the first split
1620  * criteria gives to range criteria a chance to matter, because
1621  * non-overlapping splits are equivalent in this criteria.
1622  */
1623  if (non_negative(overlap) < non_negative(context->overlap) ||
1624  (range > context->range &&
1625  non_negative(overlap) <= non_negative(context->overlap)))
1626  selectthis = true;
1627  }
1628 
1629  if (selectthis)
1630  {
1631  /* save information about selected split */
1632  context->first = false;
1633  context->ratio = ratio;
1634  context->range = range;
1635  context->overlap = overlap;
1636  context->rightLower = rightLower;
1637  context->leftUpper = leftUpper;
1638  context->dim = dimNum;
1639  POSTGIS_DEBUG(5, "split selected");
1640  }
1641  }
1642 }
1643 
1644 /*
1645  * Return increase of original BOX2DF area by new BOX2DF area insertion.
1646  */
1647 static float
1648 box_penalty(BOX2DF *original, BOX2DF *new)
1649 {
1650  float union_width,
1651  union_height;
1652 
1653  union_width = Max(original->xmax, new->xmax) -
1654  Min(original->xmin, new->xmin);
1655  union_height = Max(original->ymax, new->ymax) -
1656  Min(original->ymin, new->ymin);
1657  return union_width * union_height - (original->xmax - original->xmin) *
1658  (original->ymax - original->ymin);
1659 }
1660 
1661 /*
1662  * Compare common entries by their deltas.
1663  */
1664 static int
1665 common_entry_cmp(const void *i1, const void *i2)
1666 {
1667  float delta1 = ((const CommonEntry *) i1)->delta,
1668  delta2 = ((const CommonEntry *) i2)->delta;
1669 
1670  if (delta1 < delta2)
1671  return -1;
1672  else if (delta1 > delta2)
1673  return 1;
1674  else
1675  return 0;
1676 }
1677 
1678 /*
1679  * --------------------------------------------------------------------------
1680  * Double sorting split algorithm. This is used for both boxes and points.
1681  *
1682  * The algorithm finds split of boxes by considering splits along each axis.
1683  * Each entry is first projected as an interval on the X-axis, and different
1684  * ways to split the intervals into two groups are considered, trying to
1685  * minimize the overlap of the groups. Then the same is repeated for the
1686  * Y-axis, and the overall best split is chosen. The quality of a split is
1687  * determined by overlap along that axis and some other criteria (see
1688  * g_box_consider_split).
1689  *
1690  * After that, all the entries are divided into three groups:
1691  *
1692  * 1) Entries which should be placed to the left group
1693  * 2) Entries which should be placed to the right group
1694  * 3) "Common entries" which can be placed to any of groups without affecting
1695  * of overlap along selected axis.
1696  *
1697  * The common entries are distributed by minimizing penalty.
1698  *
1699  * For details see:
1700  * "A new double sorting-based node splitting algorithm for R-tree", A. Korotkov
1701  * http://syrcose.ispras.ru/2011/files/SYRCoSE2011_Proceedings.pdf#page=36
1702  * --------------------------------------------------------------------------
1703  */
1705 Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS)
1706 {
1707  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
1708  GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
1709  OffsetNumber i,
1710  maxoff;
1711  ConsiderSplitContext context;
1712  BOX2DF *box,
1713  *leftBox,
1714  *rightBox;
1715  int dim,
1716  commonEntriesCount;
1717  SplitInterval *intervalsLower,
1718  *intervalsUpper;
1719  CommonEntry *commonEntries;
1720  int nentries;
1721 
1722  POSTGIS_DEBUG(3, "[GIST] 'picksplit' entered");
1723 
1724  memset(&context, 0, sizeof(ConsiderSplitContext));
1725 
1726  maxoff = entryvec->n - 1;
1727  nentries = context.entriesCount = maxoff - FirstOffsetNumber + 1;
1728 
1729  /* Allocate arrays for intervals along axes */
1730  intervalsLower = (SplitInterval *) palloc(nentries * sizeof(SplitInterval));
1731  intervalsUpper = (SplitInterval *) palloc(nentries * sizeof(SplitInterval));
1732 
1733  /*
1734  * Calculate the overall minimum bounding box over all the entries.
1735  */
1736  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1737  {
1738  box = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
1739  if (i == FirstOffsetNumber)
1740  context.boundingBox = *box;
1741  else
1742  adjustBox(&context.boundingBox, box);
1743  }
1744 
1745  POSTGIS_DEBUGF(4, "boundingBox is %s", box2df_to_string(
1746  &context.boundingBox));
1747 
1748  /*
1749  * Iterate over axes for optimal split searching.
1750  */
1751  context.first = true; /* nothing selected yet */
1752  for (dim = 0; dim < 2; dim++)
1753  {
1754  float leftUpper,
1755  rightLower;
1756  int i1,
1757  i2;
1758 
1759  /* Project each entry as an interval on the selected axis. */
1760  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1761  {
1762  box = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
1763  if (dim == 0)
1764  {
1765  intervalsLower[i - FirstOffsetNumber].lower = box->xmin;
1766  intervalsLower[i - FirstOffsetNumber].upper = box->xmax;
1767  }
1768  else
1769  {
1770  intervalsLower[i - FirstOffsetNumber].lower = box->ymin;
1771  intervalsLower[i - FirstOffsetNumber].upper = box->ymax;
1772  }
1773  }
1774 
1775  /*
1776  * Make two arrays of intervals: one sorted by lower bound and another
1777  * sorted by upper bound.
1778  */
1779  memcpy(intervalsUpper, intervalsLower,
1780  sizeof(SplitInterval) * nentries);
1781  qsort(intervalsLower, nentries, sizeof(SplitInterval),
1783  qsort(intervalsUpper, nentries, sizeof(SplitInterval),
1785 
1786  /*----
1787  * The goal is to form a left and right interval, so that every entry
1788  * interval is contained by either left or right interval (or both).
1789  *
1790  * For example, with the intervals (0,1), (1,3), (2,3), (2,4):
1791  *
1792  * 0 1 2 3 4
1793  * +-+
1794  * +---+
1795  * +-+
1796  * +---+
1797  *
1798  * The left and right intervals are of the form (0,a) and (b,4).
1799  * We first consider splits where b is the lower bound of an entry.
1800  * We iterate through all entries, and for each b, calculate the
1801  * smallest possible a. Then we consider splits where a is the
1802  * uppper bound of an entry, and for each a, calculate the greatest
1803  * possible b.
1804  *
1805  * In the above example, the first loop would consider splits:
1806  * b=0: (0,1)-(0,4)
1807  * b=1: (0,1)-(1,4)
1808  * b=2: (0,3)-(2,4)
1809  *
1810  * And the second loop:
1811  * a=1: (0,1)-(1,4)
1812  * a=3: (0,3)-(2,4)
1813  * a=4: (0,4)-(2,4)
1814  */
1815 
1816  /*
1817  * Iterate over lower bound of right group, finding smallest possible
1818  * upper bound of left group.
1819  */
1820  i1 = 0;
1821  i2 = 0;
1822  rightLower = intervalsLower[i1].lower;
1823  leftUpper = intervalsUpper[i2].lower;
1824  while (true)
1825  {
1826  /*
1827  * Find next lower bound of right group.
1828  */
1829  while (i1 < nentries && (rightLower == intervalsLower[i1].lower ||
1830  isnan(intervalsLower[i1].lower)))
1831  {
1832  leftUpper = Max(leftUpper, intervalsLower[i1].upper);
1833  i1++;
1834  }
1835  if (i1 >= nentries)
1836  break;
1837  rightLower = intervalsLower[i1].lower;
1838 
1839  /*
1840  * Find count of intervals which anyway should be placed to the
1841  * left group.
1842  */
1843  while (i2 < nentries && intervalsUpper[i2].upper <= leftUpper)
1844  i2++;
1845 
1846  /*
1847  * Consider found split.
1848  */
1849  g_box_consider_split(&context, dim, rightLower, i1, leftUpper, i2);
1850  }
1851 
1852  /*
1853  * Iterate over upper bound of left group finding greates possible
1854  * lower bound of right group.
1855  */
1856  i1 = nentries - 1;
1857  i2 = nentries - 1;
1858  rightLower = intervalsLower[i1].upper;
1859  leftUpper = intervalsUpper[i2].upper;
1860  while (true)
1861  {
1862  /*
1863  * Find next upper bound of left group.
1864  */
1865  while (i2 >= 0 && (leftUpper == intervalsUpper[i2].upper ||
1866  isnan(intervalsUpper[i2].upper)))
1867  {
1868  rightLower = Min(rightLower, intervalsUpper[i2].lower);
1869  i2--;
1870  }
1871  if (i2 < 0)
1872  break;
1873  leftUpper = intervalsUpper[i2].upper;
1874 
1875  /*
1876  * Find count of intervals which anyway should be placed to the
1877  * right group.
1878  */
1879  while (i1 >= 0 && intervalsLower[i1].lower >= rightLower)
1880  i1--;
1881 
1882  /*
1883  * Consider found split.
1884  */
1885  g_box_consider_split(&context, dim,
1886  rightLower, i1 + 1, leftUpper, i2 + 1);
1887  }
1888  }
1889 
1890  /*
1891  * If we failed to find any acceptable splits, use trivial split.
1892  */
1893  if (context.first)
1894  {
1895  POSTGIS_DEBUG(4, "no acceptable splits, trivial split");
1896  fallbackSplit(entryvec, v);
1897  PG_RETURN_POINTER(v);
1898  }
1899 
1900  /*
1901  * Ok, we have now selected the split across one axis.
1902  *
1903  * While considering the splits, we already determined that there will be
1904  * enough entries in both groups to reach the desired ratio, but we did
1905  * not memorize which entries go to which group. So determine that now.
1906  */
1907 
1908  POSTGIS_DEBUGF(4, "split direction: %d", context.dim);
1909 
1910  /* Allocate vectors for results */
1911  v->spl_left = (OffsetNumber *) palloc(nentries * sizeof(OffsetNumber));
1912  v->spl_right = (OffsetNumber *) palloc(nentries * sizeof(OffsetNumber));
1913  v->spl_nleft = 0;
1914  v->spl_nright = 0;
1915 
1916  /* Allocate bounding boxes of left and right groups */
1917  leftBox = palloc0(sizeof(BOX2DF));
1918  rightBox = palloc0(sizeof(BOX2DF));
1919 
1920  /*
1921  * Allocate an array for "common entries" - entries which can be placed to
1922  * either group without affecting overlap along selected axis.
1923  */
1924  commonEntriesCount = 0;
1925  commonEntries = (CommonEntry *) palloc(nentries * sizeof(CommonEntry));
1926 
1927  /* Helper macros to place an entry in the left or right group */
1928 #define PLACE_LEFT(box, off) \
1929  do { \
1930  if (v->spl_nleft > 0) \
1931  adjustBox(leftBox, box); \
1932  else \
1933  *leftBox = *(box); \
1934  v->spl_left[v->spl_nleft++] = off; \
1935  } while(0)
1936 
1937 #define PLACE_RIGHT(box, off) \
1938  do { \
1939  if (v->spl_nright > 0) \
1940  adjustBox(rightBox, box); \
1941  else \
1942  *rightBox = *(box); \
1943  v->spl_right[v->spl_nright++] = off; \
1944  } while(0)
1945 
1946  /*
1947  * Distribute entries which can be distributed unambiguously, and collect
1948  * common entries.
1949  */
1950  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1951  {
1952  float lower,
1953  upper;
1954 
1955  /*
1956  * Get upper and lower bounds along selected axis.
1957  */
1958  box = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
1959  if (context.dim == 0)
1960  {
1961  lower = box->xmin;
1962  upper = box->xmax;
1963  }
1964  else
1965  {
1966  lower = box->ymin;
1967  upper = box->ymax;
1968  }
1969 
1970  if (upper <= context.leftUpper || isnan(upper))
1971  {
1972  /* Fits to the left group */
1973  if (lower >= context.rightLower || isnan(lower))
1974  {
1975  /* Fits also to the right group, so "common entry" */
1976  commonEntries[commonEntriesCount++].index = i;
1977  }
1978  else
1979  {
1980  /* Doesn't fit to the right group, so join to the left group */
1981  PLACE_LEFT(box, i);
1982  }
1983  }
1984  else
1985  {
1986  /*
1987  * Each entry should fit on either left or right group. Since this
1988  * entry didn't fit on the left group, it better fit in the right
1989  * group.
1990  */
1991  Assert(lower >= context.rightLower);
1992 
1993  /* Doesn't fit to the left group, so join to the right group */
1994  PLACE_RIGHT(box, i);
1995  }
1996  }
1997 
1998  POSTGIS_DEBUGF(4, "leftBox is %s", box2df_to_string(leftBox));
1999  POSTGIS_DEBUGF(4, "rightBox is %s", box2df_to_string(rightBox));
2000 
2001  /*
2002  * Distribute "common entries", if any.
2003  */
2004  if (commonEntriesCount > 0)
2005  {
2006  /*
2007  * Calculate minimum number of entries that must be placed in both
2008  * groups, to reach LIMIT_RATIO.
2009  */
2010  int m = ceil(LIMIT_RATIO * (double) nentries);
2011 
2012  /*
2013  * Calculate delta between penalties of join "common entries" to
2014  * different groups.
2015  */
2016  for (i = 0; i < commonEntriesCount; i++)
2017  {
2018  box = (BOX2DF *) DatumGetPointer(entryvec->vector[
2019  commonEntries[i].index].key);
2020  commonEntries[i].delta = Abs(box_penalty(leftBox, box) -
2021  box_penalty(rightBox, box));
2022  }
2023 
2024  /*
2025  * Sort "common entries" by calculated deltas in order to distribute
2026  * the most ambiguous entries first.
2027  */
2028  qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp);
2029 
2030  /*
2031  * Distribute "common entries" between groups.
2032  */
2033  for (i = 0; i < commonEntriesCount; i++)
2034  {
2035  box = (BOX2DF *) DatumGetPointer(entryvec->vector[
2036  commonEntries[i].index].key);
2037 
2038  /*
2039  * Check if we have to place this entry in either group to achieve
2040  * LIMIT_RATIO.
2041  */
2042  if (v->spl_nleft + (commonEntriesCount - i) <= m)
2043  PLACE_LEFT(box, commonEntries[i].index);
2044  else if (v->spl_nright + (commonEntriesCount - i) <= m)
2045  PLACE_RIGHT(box, commonEntries[i].index);
2046  else
2047  {
2048  /* Otherwise select the group by minimal penalty */
2049  if (box_penalty(leftBox, box) < box_penalty(rightBox, box))
2050  PLACE_LEFT(box, commonEntries[i].index);
2051  else
2052  PLACE_RIGHT(box, commonEntries[i].index);
2053  }
2054  }
2055  }
2056  v->spl_ldatum = PointerGetDatum(leftBox);
2057  v->spl_rdatum = PointerGetDatum(rightBox);
2058 
2059  POSTGIS_DEBUG(4, "[GIST] 'picksplit' completed");
2060 
2061  PG_RETURN_POINTER(v);
2062 }
2063 
2064 #else /* !KOROTOV_SPLIT */
2065 
2066 typedef struct
2067 {
2068  BOX2DF *key;
2069  int pos;
2070 }
2071 KBsort;
2072 
2073 static int
2074 compare_KB(const void* a, const void* b)
2075 {
2076  BOX2DF *abox = ((KBsort*)a)->key;
2077  BOX2DF *bbox = ((KBsort*)b)->key;
2078  float sa = (abox->xmax - abox->xmin) * (abox->ymax - abox->ymin);
2079  float sb = (bbox->xmax - bbox->xmin) * (bbox->ymax - bbox->ymin);
2080 
2081  if ( sa==sb ) return 0;
2082  return ( sa>sb ) ? 1 : -1;
2083 }
2084 
2091 Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS)
2092 {
2093  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
2094 
2095  GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
2096  OffsetNumber i;
2097  OffsetNumber *listL, *listR, *listB, *listT;
2098  BOX2DF *unionL, *unionR, *unionB, *unionT;
2099  int posL, posR, posB, posT;
2100  BOX2DF pageunion;
2101  BOX2DF *cur;
2102  char direction = ' ';
2103  bool allisequal = true;
2104  OffsetNumber maxoff;
2105  int nbytes;
2106 
2107  POSTGIS_DEBUG(3, "[GIST] 'picksplit' entered");
2108 
2109  posL = posR = posB = posT = 0;
2110 
2111  maxoff = entryvec->n - 1;
2112  cur = (BOX2DF*) DatumGetPointer(entryvec->vector[FirstOffsetNumber].key);
2113 
2114  memcpy((void *) &pageunion, (void *) cur, sizeof(BOX2DF));
2115 
2116  /* find MBR */
2117  for (i = OffsetNumberNext(FirstOffsetNumber); i <= maxoff; i = OffsetNumberNext(i))
2118  {
2119  cur = (BOX2DF *) DatumGetPointer(entryvec->vector[i].key);
2120 
2121  if ( allisequal == true && (
2122  pageunion.xmax != cur->xmax ||
2123  pageunion.ymax != cur->ymax ||
2124  pageunion.xmin != cur->xmin ||
2125  pageunion.ymin != cur->ymin
2126  ) )
2127  allisequal = false;
2128 
2129  if (pageunion.xmax < cur->xmax)
2130  pageunion.xmax = cur->xmax;
2131  if (pageunion.xmin > cur->xmin)
2132  pageunion.xmin = cur->xmin;
2133  if (pageunion.ymax < cur->ymax)
2134  pageunion.ymax = cur->ymax;
2135  if (pageunion.ymin > cur->ymin)
2136  pageunion.ymin = cur->ymin;
2137  }
2138 
2139  POSTGIS_DEBUGF(4, "pageunion is %s", box2df_to_string(&pageunion));
2140 
2141  nbytes = (maxoff + 2) * sizeof(OffsetNumber);
2142  listL = (OffsetNumber *) palloc(nbytes);
2143  listR = (OffsetNumber *) palloc(nbytes);
2144  unionL = (BOX2DF *) palloc(sizeof(BOX2DF));
2145  unionR = (BOX2DF *) palloc(sizeof(BOX2DF));
2146 
2147  if (allisequal)
2148  {
2149  POSTGIS_DEBUG(4, " AllIsEqual!");
2150 
2151  cur = (BOX2DF*) DatumGetPointer(entryvec->vector[OffsetNumberNext(FirstOffsetNumber)].key);
2152 
2153 
2154  if (memcmp((void *) cur, (void *) &pageunion, sizeof(BOX2DF)) == 0)
2155  {
2156  v->spl_left = listL;
2157  v->spl_right = listR;
2158  v->spl_nleft = v->spl_nright = 0;
2159  memcpy((void *) unionL, (void *) &pageunion, sizeof(BOX2DF));
2160  memcpy((void *) unionR, (void *) &pageunion, sizeof(BOX2DF));
2161 
2162  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
2163  {
2164  if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
2165  {
2166  v->spl_left[v->spl_nleft] = i;
2167  v->spl_nleft++;
2168  }
2169  else
2170  {
2171  v->spl_right[v->spl_nright] = i;
2172  v->spl_nright++;
2173  }
2174  }
2175  v->spl_ldatum = PointerGetDatum(unionL);
2176  v->spl_rdatum = PointerGetDatum(unionR);
2177 
2178  PG_RETURN_POINTER(v);
2179  }
2180  }
2181 
2182  listB = (OffsetNumber *) palloc(nbytes);
2183  listT = (OffsetNumber *) palloc(nbytes);
2184  unionB = (BOX2DF *) palloc(sizeof(BOX2DF));
2185  unionT = (BOX2DF *) palloc(sizeof(BOX2DF));
2186 
2187 #define ADDLIST( list, unionD, pos, num ) do { \
2188  if ( pos ) { \
2189  if ( unionD->xmax < cur->xmax ) unionD->xmax = cur->xmax; \
2190  if ( unionD->xmin > cur->xmin ) unionD->xmin = cur->xmin; \
2191  if ( unionD->ymax < cur->ymax ) unionD->ymax = cur->ymax; \
2192  if ( unionD->ymin > cur->ymin ) unionD->ymin = cur->ymin; \
2193  } else { \
2194  memcpy( (void*)unionD, (void*) cur, sizeof( BOX2DF ) ); \
2195  } \
2196  list[pos] = num; \
2197  (pos)++; \
2198 } while(0)
2199 
2200  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
2201  {
2202  cur = (BOX2DF*) DatumGetPointer(entryvec->vector[i].key);
2203 
2204  if (cur->xmin - pageunion.xmin < pageunion.xmax - cur->xmax)
2205  ADDLIST(listL, unionL, posL,i);
2206  else
2207  ADDLIST(listR, unionR, posR,i);
2208  if (cur->ymin - pageunion.ymin < pageunion.ymax - cur->ymax)
2209  ADDLIST(listB, unionB, posB,i);
2210  else
2211  ADDLIST(listT, unionT, posT,i);
2212  }
2213 
2214  POSTGIS_DEBUGF(4, "unionL is %s", box2df_to_string(unionL));
2215  POSTGIS_DEBUGF(4, "unionR is %s", box2df_to_string(unionR));
2216  POSTGIS_DEBUGF(4, "unionT is %s", box2df_to_string(unionT));
2217  POSTGIS_DEBUGF(4, "unionB is %s", box2df_to_string(unionB));
2218 
2219  /* bad disposition, sort by ascending and resplit */
2220  if ( (posR==0 || posL==0) && (posT==0 || posB==0) )
2221  {
2222  KBsort *arr = (KBsort*)palloc( sizeof(KBsort) * maxoff );
2223  posL = posR = posB = posT = 0;
2224  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
2225  {
2226  arr[i-1].key = (BOX2DF*) DatumGetPointer(entryvec->vector[i].key);
2227  arr[i-1].pos = i;
2228  }
2229  qsort( arr, maxoff, sizeof(KBsort), compare_KB );
2230  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
2231  {
2232  cur = arr[i-1].key;
2233  if (cur->xmin - pageunion.xmin < pageunion.xmax - cur->xmax)
2234  ADDLIST(listL, unionL, posL,arr[i-1].pos);
2235  else if ( cur->xmin - pageunion.xmin == pageunion.xmax - cur->xmax )
2236  {
2237  if ( posL>posR )
2238  ADDLIST(listR, unionR, posR,arr[i-1].pos);
2239  else
2240  ADDLIST(listL, unionL, posL,arr[i-1].pos);
2241  }
2242  else
2243  ADDLIST(listR, unionR, posR,arr[i-1].pos);
2244 
2245  if (cur->ymin - pageunion.ymin < pageunion.ymax - cur->ymax)
2246  ADDLIST(listB, unionB, posB,arr[i-1].pos);
2247  else if ( cur->ymin - pageunion.ymin == pageunion.ymax - cur->ymax )
2248  {
2249  if ( posB>posT )
2250  ADDLIST(listT, unionT, posT,arr[i-1].pos);
2251  else
2252  ADDLIST(listB, unionB, posB,arr[i-1].pos);
2253  }
2254  else
2255  ADDLIST(listT, unionT, posT,arr[i-1].pos);
2256  }
2257  pfree(arr);
2258  }
2259 
2260  /* which split more optimal? */
2261  if (Max(posL, posR) < Max(posB, posT))
2262  direction = 'x';
2263  else if (Max(posL, posR) > Max(posB, posT))
2264  direction = 'y';
2265  else
2266  {
2267  float sizeLR, sizeBT;
2268  BOX2DF interLR, interBT;
2269 
2270  if ( box2df_intersection(unionL, unionR, &interLR) == FALSE )
2271  sizeLR = 0.0;
2272  else
2273  sizeLR = box2df_size(&interLR);
2274 
2275  if ( box2df_intersection(unionB, unionT, &interBT) == FALSE )
2276  sizeBT = 0.0;
2277  else
2278  sizeBT = box2df_size(&interBT);
2279 
2280  if (sizeLR < sizeBT)
2281  direction = 'x';
2282  else
2283  direction = 'y';
2284  }
2285 
2286  POSTGIS_DEBUGF(4, "split direction '%c'", direction);
2287 
2288  if (direction == 'x')
2289  {
2290  pfree(unionB);
2291  pfree(listB);
2292  pfree(unionT);
2293  pfree(listT);
2294 
2295  v->spl_left = listL;
2296  v->spl_right = listR;
2297  v->spl_nleft = posL;
2298  v->spl_nright = posR;
2299  v->spl_ldatum = PointerGetDatum(unionL);
2300  v->spl_rdatum = PointerGetDatum(unionR);
2301  }
2302  else
2303  {
2304  pfree(unionR);
2305  pfree(listR);
2306  pfree(unionL);
2307  pfree(listL);
2308 
2309  v->spl_left = listB;
2310  v->spl_right = listT;
2311  v->spl_nleft = posB;
2312  v->spl_nright = posT;
2313  v->spl_ldatum = PointerGetDatum(unionB);
2314  v->spl_rdatum = PointerGetDatum(unionT);
2315  }
2316 
2317  POSTGIS_DEBUG(4, "[GIST] 'picksplit' completed");
2318 
2319  PG_RETURN_POINTER(v);
2320 }
2321 
2322 #endif
2323 
2324 
2325 /*
2326 ** The BOX32DF key must be defined as a PostgreSQL type, even though it is only
2327 ** ever used internally. These no-op stubs are used to bind the type.
2328 */
2330 Datum box2df_in(PG_FUNCTION_ARGS)
2331 {
2332  ereport(ERROR,(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2333  errmsg("function box2df_in not implemented")));
2334  PG_RETURN_POINTER(NULL);
2335 }
2336 
2338 Datum box2df_out(PG_FUNCTION_ARGS)
2339 {
2340  BOX2DF *box = (BOX2DF *) PG_GETARG_POINTER(0);
2341  char *result = box2df_to_string(box);
2342  PG_RETURN_CSTRING(result);
2343 }
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *box)
Read the bounding box off a serialization and calculate one if it is not already there.
Definition: g_serialized.c:398
static float non_negative(float val)
Datum gserialized_above_2d(PG_FUNCTION_ARGS)
Datum gserialized_overabove_2d(PG_FUNCTION_ARGS)
static void box2df_validate(BOX2DF *b)
uint8_t data[1]
Definition: liblwgeom.h:383
#define LIMIT_RATIO
Datum gserialized_gist_penalty_2d(PG_FUNCTION_ARGS)
static BOX2DF * box2df_copy(BOX2DF *b)
static bool box2df_overlaps(const BOX2DF *a, const BOX2DF *b)
static double box2df_distance_node_centroid(const BOX2DF *node, const BOX2DF *query)
Calculate the The node_box_edge->query_centroid distance between the boxes.
Datum gserialized_gist_compress_2d(PG_FUNCTION_ARGS)
static bool box2df_left(const BOX2DF *a, const BOX2DF *b)
double xmax
Definition: liblwgeom.h:292
static void fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
#define LW_SUCCESS
Definition: liblwgeom.h:79
static bool box2df_below(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_left_2d(PG_FUNCTION_ARGS)
static bool gserialized_gist_consistent_leaf_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
static void g_box_consider_split(ConsiderSplitContext *context, int dimNum, float rightLower, int minLeftCount, float leftUpper, int maxLeftCount)
static int interval_cmp_lower(const void *i1, const void *i2)
static void adjustBox(BOX2DF *b, BOX2DF *addon)
static bool box2df_right(const BOX2DF *a, const BOX2DF *b)
float next_float_down(double d)
Definition: lwgeom_api.c:156
static bool box2df_above(const BOX2DF *a, const BOX2DF *b)
static float box2df_size(const BOX2DF *a)
bool(* box2df_predicate)(const BOX2DF *a, const BOX2DF *b)
static bool box2df_within(const BOX2DF *a, const BOX2DF *b)
#define LW_FAILURE
Definition: liblwgeom.h:78
float next_float_up(double d)
Definition: lwgeom_api.c:172
static int box2df_from_gbox_p(GBOX *box, BOX2DF *a)
double ymin
Definition: liblwgeom.h:293
Datum gserialized_overlaps_2d(PG_FUNCTION_ARGS)
Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS)
Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS)
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:51
double xmin
Definition: liblwgeom.h:291
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.
#define FLAGS_GET_BBOX(flags)
Definition: liblwgeom.h:141
static char * box2df_to_string(const BOX2DF *a)
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
Datum gserialized_below_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_same_2d(PG_FUNCTION_ARGS)
static void box2df_merge(BOX2DF *b_union, BOX2DF *b_new)
Datum gserialized_within_2d(PG_FUNCTION_ARGS)
int gserialized_datum_get_box2df_p(Datum gsdatum, BOX2DF *box2df)
Peak into a GSERIALIZED datum to find the bounding box.
Datum gserialized_overleft_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS)
Datum box2df_out(PG_FUNCTION_ARGS)
PG_FUNCTION_INFO_V1(gserialized_distance_centroid_2d)
double ymax
Definition: liblwgeom.h:294
Datum gserialized_same_2d(PG_FUNCTION_ARGS)
static int interval_cmp_upper(const void *i1, const void *i2)
static double pt_distance(double ax, double ay, double bx, double by)
Datum distance(PG_FUNCTION_ARGS)
static bool box2df_overright(const BOX2DF *a, const BOX2DF *b)
static int common_entry_cmp(const void *i1, const void *i2)
Datum gserialized_gist_picksplit_2d(PG_FUNCTION_ARGS)
static bool box2df_overabove(const BOX2DF *a, const BOX2DF *b)
static bool box2df_overbelow(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_gist_decompress_2d(PG_FUNCTION_ARGS)
static bool box2df_overleft(const BOX2DF *a, const BOX2DF *b)
#define FALSE
Definition: dbfopen.c:168
Datum gserialized_overright_2d(PG_FUNCTION_ARGS)
bool box2df_contains(const BOX2DF *a, const BOX2DF *b)
#define PLACE_RIGHT(box, off)
static float box2df_union_size(const BOX2DF *a, const BOX2DF *b)
static int gserialized_datum_predicate_2d(Datum gs1, Datum gs2, box2df_predicate predicate)
Support function.
static bool box2df_equals(const BOX2DF *a, const BOX2DF *b)
Datum gserialized_right_2d(PG_FUNCTION_ARGS)
Datum box2df_in(PG_FUNCTION_ARGS)
#define TRUE
Definition: dbfopen.c:169
Datum gserialized_contains_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_union_2d(PG_FUNCTION_ARGS)
Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
#define PLACE_LEFT(box, off)
static bool gserialized_gist_consistent_internal_2d(BOX2DF *key, BOX2DF *query, StrategyNumber strategy)
static float box_penalty(BOX2DF *original, BOX2DF *new)
Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS)
This library is the generic geometry handling section of PostGIS.
uint8_t flags
Definition: liblwgeom.h:382