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