PostGIS  2.3.7dev-r@@SVN_REVISION@@
gserialized_gist_nd.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 */
38 
39 #include "postgres.h"
40 #include "access/gist.h" /* For GiST */
41 #include "access/itup.h"
42 #include "access/skey.h"
43 
44 #include "../postgis_config.h"
45 
46 /*#define POSTGIS_DEBUG_LEVEL 4*/
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 #include "geography.h"
52 
53 #include <assert.h>
54 
55 #if POSTGIS_PGSQL_VERSION >= 100
56 #include <float.h>
57 #endif
58 
59 /* Fall back to older finite() if necessary */
60 #ifndef HAVE_ISFINITE
61 # ifdef HAVE_GNU_ISFINITE
62 # define _GNU_SOURCE
63 # else
64 # define isfinite finite
65 # endif
66 #endif
67 
68 /*
69 ** When is a node split not so good? If more than 90% of the entries
70 ** end up in one of the children.
71 */
72 #define LIMIT_RATIO 0.1
73 
74 /*
75 ** For debugging
76 */
77 #if POSTGIS_DEBUG_LEVEL > 0
78 static int geog_counter_leaf = 0;
79 static int geog_counter_internal = 0;
80 #endif
81 
82 /*
83 ** ND Index key type stub prototypes
84 */
85 Datum gidx_out(PG_FUNCTION_ARGS);
86 Datum gidx_in(PG_FUNCTION_ARGS);
87 
88 /*
89 ** ND GiST prototypes
90 */
91 Datum gserialized_gist_consistent(PG_FUNCTION_ARGS);
92 Datum gserialized_gist_compress(PG_FUNCTION_ARGS);
93 Datum gserialized_gist_decompress(PG_FUNCTION_ARGS);
94 Datum gserialized_gist_penalty(PG_FUNCTION_ARGS);
95 Datum gserialized_gist_picksplit(PG_FUNCTION_ARGS);
96 Datum gserialized_gist_union(PG_FUNCTION_ARGS);
97 Datum gserialized_gist_same(PG_FUNCTION_ARGS);
98 Datum gserialized_gist_distance(PG_FUNCTION_ARGS);
99 Datum gserialized_gist_geog_distance(PG_FUNCTION_ARGS);
100 
101 /*
102 ** ND Operator prototypes
103 */
104 Datum gserialized_overlaps(PG_FUNCTION_ARGS);
105 Datum gserialized_contains(PG_FUNCTION_ARGS);
106 #if POSTGIS_PGSQL_VERSION > 94
107 Datum gserialized_gidx_geom_contains(PG_FUNCTION_ARGS);
108 Datum gserialized_gidx_gidx_contains(PG_FUNCTION_ARGS);
109 #endif
110 Datum gserialized_within(PG_FUNCTION_ARGS);
111 #if POSTGIS_PGSQL_VERSION > 94
112 Datum gserialized_gidx_geom_within(PG_FUNCTION_ARGS);
113 Datum gserialized_gidx_gidx_within(PG_FUNCTION_ARGS);
114 #endif
115 Datum gserialized_distance_nd(PG_FUNCTION_ARGS);
116 #if POSTGIS_PGSQL_VERSION > 94
117 Datum gserialized_gidx_geom_same(PG_FUNCTION_ARGS);
118 Datum gserialized_gidx_gidx_same(PG_FUNCTION_ARGS);
119 #endif
120 
121 /*
122 ** GIDX true/false test function type
123 */
124 typedef bool (*gidx_predicate)(GIDX *a, GIDX *b);
125 
126 
127 /* Allocate a new copy of GIDX */
128 static GIDX* gidx_copy(GIDX *b)
129 {
130  GIDX *c = (GIDX*)palloc(VARSIZE(b));
131  POSTGIS_DEBUGF(5, "copied gidx (%p) to gidx (%p)", b, c);
132  memcpy((void*)c, (void*)b, VARSIZE(b));
133  return c;
134 }
135 
136 
137 /* Ensure all minimums are below maximums. */
138 static inline void gidx_validate(GIDX *b)
139 {
140  int i;
141  Assert(b);
142  POSTGIS_DEBUGF(5,"validating gidx (%s)", gidx_to_string(b));
143  for ( i = 0; i < GIDX_NDIMS(b); i++ )
144  {
145  if ( GIDX_GET_MIN(b,i) > GIDX_GET_MAX(b,i) )
146  {
147  float tmp;
148  tmp = GIDX_GET_MIN(b,i);
149  GIDX_SET_MIN(b,i,GIDX_GET_MAX(b,i));
150  GIDX_SET_MAX(b,i,tmp);
151  }
152  }
153  return;
154 }
155 
156 /* An "unknown" GIDX is used to represent the bounds of an EMPTY
157  geometry or other-wise unindexable geometry (like one with NaN
158  or Inf bounds) */
159 static inline bool gidx_is_unknown(const GIDX *a)
160 {
161  size_t size = VARSIZE(a) - VARHDRSZ;
162  /* "unknown" gidx objects have a too-small size of one float */
163  if ( size <= 0.0 )
164  return TRUE;
165  return FALSE;
166 }
167 
168 static inline void gidx_set_unknown(GIDX *a)
169 {
170  SET_VARSIZE(a, VARHDRSZ);
171 }
172 
173 /* Enlarge b_union to contain b_new. If b_new contains more
174  dimensions than b_union, expand b_union to contain those dimensions. */
175 static void gidx_merge(GIDX **b_union, GIDX *b_new)
176 {
177  int i, dims_union, dims_new;
178  Assert(b_union);
179  Assert(*b_union);
180  Assert(b_new);
181 
182  /* Can't merge an unknown into any thing */
183  if( gidx_is_unknown(b_new) )
184  return;
185 
186  /* Merge of unknown and known is known */
187  if( gidx_is_unknown(*b_union) )
188  {
189  *b_union = b_new;
190  return;
191  }
192 
193  dims_union = GIDX_NDIMS(*b_union);
194  dims_new = GIDX_NDIMS(b_new);
195 
196  POSTGIS_DEBUGF(4, "merging gidx (%s) into gidx (%s)", gidx_to_string(b_new), gidx_to_string(*b_union));
197 
198  if ( dims_new > dims_union )
199  {
200  POSTGIS_DEBUGF(5, "reallocating b_union from %d dims to %d dims", dims_union, dims_new);
201  *b_union = (GIDX*)repalloc(*b_union, GIDX_SIZE(dims_new));
202  SET_VARSIZE(*b_union, VARSIZE(b_new));
203  dims_union = dims_new;
204  }
205 
206  for ( i = 0; i < dims_new; i++ )
207  {
208  /* Adjust minimums */
209  GIDX_SET_MIN(*b_union, i, Min(GIDX_GET_MIN(*b_union,i),GIDX_GET_MIN(b_new,i)));
210  /* Adjust maximums */
211  GIDX_SET_MAX(*b_union, i, Max(GIDX_GET_MAX(*b_union,i),GIDX_GET_MAX(b_new,i)));
212  }
213 
214  POSTGIS_DEBUGF(5, "merge complete (%s)", gidx_to_string(*b_union));
215  return;
216 }
217 
218 /* Calculate the volume (in n-d units) of the GIDX */
219 static float gidx_volume(GIDX *a)
220 {
221  float result;
222  int i;
223  if ( a == NULL || gidx_is_unknown(a) )
224  {
225  /* elog(ERROR, "gidx_volume received a null argument"); */
226  return 0.0;
227  }
228  result = GIDX_GET_MAX(a,0) - GIDX_GET_MIN(a,0);
229  for ( i = 1; i < GIDX_NDIMS(a); i++ )
230  result *= (GIDX_GET_MAX(a,i) - GIDX_GET_MIN(a,i));
231  POSTGIS_DEBUGF(5, "calculated volume of %s as %.8g", gidx_to_string(a), result);
232  return result;
233 }
234 
235 /* Ensure the first argument has the higher dimensionality. */
236 static void gidx_dimensionality_check(GIDX **a, GIDX **b)
237 {
238  if ( GIDX_NDIMS(*a) < GIDX_NDIMS(*b) )
239  {
240  GIDX *tmp = *b;
241  *b = *a;
242  *a = tmp;
243  }
244 }
245 
246 /* Calculate the volume of the union of the boxes. Avoids creating an intermediate box. */
247 static float gidx_union_volume(GIDX *a, GIDX *b)
248 {
249  float result;
250  int i;
251  int ndims_a, ndims_b;
252 
253  POSTGIS_DEBUG(5,"entered function");
254 
255  if ( a == NULL && b == NULL )
256  {
257  elog(ERROR, "gidx_union_volume received two null arguments");
258  return 0.0;
259  }
260 
261  if ( a == NULL || gidx_is_unknown(a) )
262  return gidx_volume(b);
263 
264  if ( b == NULL || gidx_is_unknown(b) )
265  return gidx_volume(a);
266 
267  if ( gidx_is_unknown(a) && gidx_is_unknown(b) )
268  {
269  return 0.0;
270  }
271 
272  /* Ensure 'a' has the most dimensions. */
274 
275  ndims_a = GIDX_NDIMS(a);
276  ndims_b = GIDX_NDIMS(b);
277 
278  /* Initialize with maximal length of first dimension. */
279  result = Max(GIDX_GET_MAX(a,0),GIDX_GET_MAX(b,0)) - Min(GIDX_GET_MIN(a,0),GIDX_GET_MIN(b,0));
280 
281  /* Multiply by maximal length of remaining dimensions. */
282  for ( i = 1; i < ndims_b; i++ )
283  {
284  result *= (Max(GIDX_GET_MAX(a,i),GIDX_GET_MAX(b,i)) - Min(GIDX_GET_MIN(a,i),GIDX_GET_MIN(b,i)));
285  }
286 
287  /* Add in dimensions of higher dimensional box. */
288  for ( i = ndims_b; i < ndims_a; i++ )
289  {
290  result *= (GIDX_GET_MAX(a,i) - GIDX_GET_MIN(a,i));
291  }
292 
293  POSTGIS_DEBUGF(5, "volume( %s union %s ) = %.12g", gidx_to_string(a), gidx_to_string(b), result);
294 
295  return result;
296 }
297 
298 /* Calculate the volume of the intersection of the boxes. */
299 static float gidx_inter_volume(GIDX *a, GIDX *b)
300 {
301  int i;
302  float result;
303 
304  POSTGIS_DEBUG(5,"entered function");
305 
306  if ( a == NULL || b == NULL )
307  {
308  elog(ERROR, "gidx_inter_volume received a null argument");
309  return 0.0;
310  }
311 
312  if ( gidx_is_unknown(a) || gidx_is_unknown(b) )
313  {
314  return 0.0;
315  }
316 
317  /* Ensure 'a' has the most dimensions. */
319 
320  /* Initialize with minimal length of first dimension. */
321  result = Min(GIDX_GET_MAX(a,0),GIDX_GET_MAX(b,0)) - Max(GIDX_GET_MIN(a,0),GIDX_GET_MIN(b,0));
322 
323  /* If they are disjoint (max < min) then return zero. */
324  if ( result < 0.0 ) return 0.0;
325 
326  /* Continue for remaining dimensions. */
327  for ( i = 1; i < GIDX_NDIMS(b); i++ )
328  {
329  float width = Min(GIDX_GET_MAX(a,i),GIDX_GET_MAX(b,i)) - Max(GIDX_GET_MIN(a,i),GIDX_GET_MIN(b,i));
330  if ( width < 0.0 ) return 0.0;
331  /* Multiply by minimal length of remaining dimensions. */
332  result *= width;
333  }
334  POSTGIS_DEBUGF(5, "volume( %s intersection %s ) = %.12g", gidx_to_string(a), gidx_to_string(b), result);
335  return result;
336 }
337 
338 /*
339 ** Overlapping GIDX box test.
340 **
341 ** Box(A) Overlaps Box(B) IFF for every dimension d:
342 ** min(A,d) <= max(B,d) && max(A,d) => min(B,d)
343 **
344 ** Any missing dimension is assumed by convention to
345 ** overlap whatever finite range available on the
346 ** other operand. See
347 ** http://lists.osgeo.org/pipermail/postgis-devel/2015-February/024757.html
348 **
349 ** Empty boxes never overlap.
350 */
351 static bool gidx_overlaps(GIDX *a, GIDX *b)
352 {
353  int i;
354  int ndims_b;
355  POSTGIS_DEBUG(5, "entered function");
356 
357  if ( (a == NULL) || (b == NULL) ) return FALSE;
358 
359  if ( gidx_is_unknown(a) || gidx_is_unknown(b) )
360  return FALSE;
361 
362  /* Ensure 'a' has the most dimensions. */
364 
365  ndims_b = GIDX_NDIMS(b);
366 
367  /* compare only up to dimensions of (b), missing dimensions always overlap */
368  for ( i = 0; i < ndims_b; i++ )
369  {
370  if ( GIDX_GET_MIN(a,i) > GIDX_GET_MAX(b,i) )
371  return FALSE;
372  if ( GIDX_GET_MIN(b,i) > GIDX_GET_MAX(a,i) )
373  return FALSE;
374  }
375 
376  return TRUE;
377 }
378 
379 /*
380 ** Containment GIDX test.
381 **
382 ** Box(A) CONTAINS Box(B) IFF (pt(A)LL < pt(B)LL) && (pt(A)UR > pt(B)UR)
383 */
384 bool gidx_contains(GIDX *a, GIDX *b)
385 {
386  int i, dims_a, dims_b;
387 
388  POSTGIS_DEBUG(5, "entered function");
389 
390  if ( (a == NULL) || (b == NULL) ) return FALSE;
391 
392  if ( gidx_is_unknown(a) || gidx_is_unknown(b) )
393  return FALSE;
394 
395  dims_a = GIDX_NDIMS(a);
396  dims_b = GIDX_NDIMS(b);
397 
398  if ( dims_a < dims_b )
399  {
400  /*
401  ** If (b) is of higher dimensionality than (a) it can only be contained
402  ** if those higher dimensions are zeroes.
403  */
404  for (i = dims_a; i < dims_b; i++)
405  {
406  if ( GIDX_GET_MIN(b,i) != 0 )
407  return FALSE;
408  if ( GIDX_GET_MAX(b,i) != 0 )
409  return FALSE;
410  }
411  }
412 
413  /* Excess dimensions of (a), don't matter, it just has to contain (b) in (b)'s dimensions */
414  for (i = 0; i < Min(dims_a, dims_b); i++)
415  {
416  if ( GIDX_GET_MIN(a,i) > GIDX_GET_MIN(b,i) )
417  return FALSE;
418  if ( GIDX_GET_MAX(a,i) < GIDX_GET_MAX(b,i) )
419  return FALSE;
420  }
421 
422  return TRUE;
423 }
424 
425 /*
426 ** Equality GIDX test.
427 **
428 ** Box(A) EQUALS Box(B) IFF (pt(A)LL == pt(B)LL) && (pt(A)UR == pt(B)UR)
429 */
430 static bool gidx_equals(GIDX *a, GIDX *b)
431 {
432  int i;
433 
434  POSTGIS_DEBUG(5, "entered function");
435 
436  if ( (a == NULL) && (b == NULL) ) return TRUE;
437  if ( (a == NULL) || (b == NULL) ) return FALSE;
438 
439  if ( gidx_is_unknown(a) && gidx_is_unknown(b) )
440  return TRUE;
441 
442  if ( gidx_is_unknown(a) || gidx_is_unknown(b) )
443  return FALSE;
444 
445  /* Ensure 'a' has the most dimensions. */
447 
448  /* For all shared dimensions min(a) == min(b), max(a) == max(b) */
449  for (i = 0; i < GIDX_NDIMS(b); i++)
450  {
451  if ( GIDX_GET_MIN(a,i) != GIDX_GET_MIN(b,i) )
452  return FALSE;
453  if ( GIDX_GET_MAX(a,i) != GIDX_GET_MAX(b,i) )
454  return FALSE;
455  }
456  /* For all unshared dimensions min(a) == 0.0, max(a) == 0.0 */
457  for (i = GIDX_NDIMS(b); i < GIDX_NDIMS(a); i++)
458  {
459  if ( GIDX_GET_MIN(a,i) != 0.0 )
460  return FALSE;
461  if ( GIDX_GET_MAX(a,i) != 0.0 )
462  return FALSE;
463  }
464  return TRUE;
465 }
466 
471 static int
472 gserialized_datum_predicate(Datum gs1, Datum gs2, gidx_predicate predicate)
473 {
474  /* Put aside some stack memory and use it for GIDX pointers. */
475  char boxmem1[GIDX_MAX_SIZE];
476  char boxmem2[GIDX_MAX_SIZE];
477  GIDX *gidx1 = (GIDX*)boxmem1;
478  GIDX *gidx2 = (GIDX*)boxmem2;
479 
480  POSTGIS_DEBUG(3, "entered function");
481 
482  /* Must be able to build box for each arguement (ie, not empty geometry)
483  and predicate function to return true. */
484  if ( (gserialized_datum_get_gidx_p(gs1, gidx1) == LW_SUCCESS) &&
485  (gserialized_datum_get_gidx_p(gs2, gidx2) == LW_SUCCESS) &&
486  predicate(gidx1, gidx2) )
487  {
488  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(gidx1), gidx_to_string(gidx2));
489  return LW_TRUE;
490  }
491  return LW_FALSE;
492 }
493 
494 #if POSTGIS_PGSQL_VERSION > 94
495 static int
496 gserialized_datum_predicate_gidx_geom(GIDX *gidx1, Datum gs2, gidx_predicate predicate)
497 {
498  /* Put aside some stack memory and use it for GIDX pointers. */
499  char boxmem2[GIDX_MAX_SIZE];
500  GIDX *gidx2 = (GIDX*)boxmem2;
501 
502  POSTGIS_DEBUG(3, "entered function");
503 
504  /* Must be able to build box for gs2 arguement (ie, not empty geometry)
505  and predicate function to return true. */
506  if ( (gserialized_datum_get_gidx_p(gs2, gidx2) == LW_SUCCESS) &&
507  predicate(gidx1, gidx2) )
508  {
509  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(gidx1), gidx_to_string(gidx2));
510  return LW_TRUE;
511  }
512  return LW_FALSE;
513 }
514 
515 static int
516 gserialized_datum_predicate_geom_gidx(Datum gs1, GIDX *gidx2, gidx_predicate predicate)
517 {
518  /* Put aside some stack memory and use it for GIDX pointers. */
519  char boxmem2[GIDX_MAX_SIZE];
520  GIDX *gidx1 = (GIDX*)boxmem2;
521 
522  POSTGIS_DEBUG(3, "entered function");
523 
524  /* Must be able to build box for gs2 arguement (ie, not empty geometry)
525  and predicate function to return true. */
526  if ( (gserialized_datum_get_gidx_p(gs1, gidx1) == LW_SUCCESS) &&
527  predicate(gidx1, gidx2) )
528  {
529  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(gidx1), gidx_to_string(gidx2));
530  return LW_TRUE;
531  }
532  return LW_FALSE;
533 }
534 #endif
535 
536 
540 #if POSTGIS_PGSQL_VERSION < 95
541 static double gidx_distance_leaf_centroid(const GIDX *a, const GIDX *b)
542 {
543  int ndims, i;
544  double sum = 0;
545 
546  /* Base computation on least available dimensions */
547  ndims = Min(GIDX_NDIMS(b), GIDX_NDIMS(a));
548  for ( i = 0; i < ndims; ++i )
549  {
550  double ca, cb, d;
551  double amin = GIDX_GET_MIN(a,i);
552  double amax = GIDX_GET_MAX(a,i);
553  double bmin = GIDX_GET_MIN(b,i);
554  double bmax = GIDX_GET_MAX(b,i);
555  ca = amin + ( ( amax - amin ) / 2.0 );
556  cb = bmin + ( ( bmax - bmin ) / 2.0 );
557  d = ca - cb;
558  if ( ! isfinite(d) )
559  {
560  /* Can happen if a dimension was padded with FLT_MAX,
561  * effectively meaning "infinite range". In that case
562  * we take that dimension as adding 0 to the total
563  * distance.
564  */
565  continue;
566  }
567  sum += d * d;
568 /*
569  POSTGIS_DEBUGF(3, " centroid of A for dimension %d is %g", i, ca);
570  POSTGIS_DEBUGF(3, " centroid of B for dimension %d is %g", i, cb);
571  POSTGIS_DEBUGF(3, " distance on dimension %d is %g, squared as %g, grows sum to %g", i, d, d*d, sum);
572 */
573  }
574  return sqrt(sum);
575 }
576 #endif
577 
581 static double gidx_distance(const GIDX *a, const GIDX *b, int m_is_time)
582 {
583  int ndims, i;
584  double sum = 0;
585 
586  /* Base computation on least available dimensions */
587  ndims = Min(GIDX_NDIMS(b), GIDX_NDIMS(a));
588  for ( i = 0; i < ndims; ++i )
589  {
590  double d;
591  double amin = GIDX_GET_MIN(a,i);
592  double amax = GIDX_GET_MAX(a,i);
593  double bmin = GIDX_GET_MIN(b,i);
594  double bmax = GIDX_GET_MAX(b,i);
595  POSTGIS_DEBUGF(3, "A %g - %g", amin, amax);
596  POSTGIS_DEBUGF(3, "B %g - %g", bmin, bmax);
597 
598  if ( ( amin <= bmax && amax >= bmin ) )
599  {
600  /* overlaps */
601  d = 0;
602  }
603  else if ( i == 4 && m_is_time )
604  {
605  return FLT_MAX;
606  }
607  else if ( bmax < amin )
608  {
609  /* is "left" */
610  d = amin - bmax;
611  }
612  else
613  {
614  /* is "right" */
615  assert( bmin > amax );
616  d = bmin - amax;
617  }
618  if ( ! isfinite(d) )
619  {
620  /* Can happen if coordinates are corrupted/NaN */
621  continue;
622  }
623  sum += d * d;
624  POSTGIS_DEBUGF(3, "dist %g, squared %g, grows sum to %g", d, d*d, sum);
625  }
626  return sqrt(sum);
627 }
628 
629 #if POSTGIS_PGSQL_VERSION < 95
630 static double gidx_distance_node_centroid(const GIDX *node, const GIDX *query)
631 {
632  int i;
633  double sum = 0;
634 
635  /* Base computation on least available dimensions */
636  int ndims = Min(GIDX_NDIMS(node), GIDX_NDIMS(query));
637 
638  for ( i = 0; i < ndims; ++i )
639  {
640  double d;
641  double amin = GIDX_GET_MIN(query,i);
642  double amax = GIDX_GET_MAX(query,i);
643  double bmin = GIDX_GET_MIN(node,i);
644  double bmax = GIDX_GET_MAX(node,i);
645  double ca = amin + ( ( amax - amin ) / 2.0 );
646 
647  if ( ( ca <= bmax && ca >= bmin ) )
648  {
649  /* overlaps */
650  d = 0;
651  }
652  else if ( bmax < ca )
653  {
654  /* is "left" */
655  d = ca - bmax;
656  }
657  else
658  {
659  /* is "right" */
660  assert( bmin > ca );
661  d = bmin - ca;
662  }
663  if ( ! isfinite(d) )
664  {
665  /* Can happen if coordinates are corrupted/NaN */
666  continue;
667  }
668  sum += d * d;
669  POSTGIS_DEBUGF(3, "dist %g, squared %g, grows sum to %g", d, d*d, sum);
670  }
671  return sqrt(sum);
672 }
673 #else /* POSTGIS_PGSQL_VERSION >= 95 */
674 static double gidx_distance_m(const GIDX *a, const GIDX *b)
675 {
676  int mdim_a, mdim_b;
677  double d, amin, amax, bmin, bmax;
678 
679  /* Base computation on least available dimensions */
680  mdim_a = GIDX_NDIMS(a) - 1;
681  mdim_b = GIDX_NDIMS(b) - 1;
682 
683  amin = GIDX_GET_MIN(a,mdim_a);
684  amax = GIDX_GET_MAX(a,mdim_a);
685  bmin = GIDX_GET_MIN(b,mdim_b);
686  bmax = GIDX_GET_MAX(b,mdim_b);
687 
688  if ( ( amin <= bmax && amax >= bmin ) )
689  {
690  /* overlaps */
691  d = 0;
692  }
693  else if ( bmax < amin )
694  {
695  /* is "left" */
696  d = amin - bmax;
697  }
698  else
699  {
700  /* is "right" */
701  assert( bmin > amax );
702  d = bmin - amax;
703  }
704 
705  return d;
706 }
707 #endif /* POSTGIS_PGSQL_VERSION >= 96 */
708 
714 {
715  char boxmem[GIDX_MAX_SIZE];
716  GIDX *gidx = (GIDX*)boxmem;
717  float fdistance = (float)distance;
718 
719  /* Get our bounding box out of the geography, return right away if
720  input is an EMPTY geometry. */
721  if ( gserialized_get_gidx_p(g, gidx) == LW_FAILURE )
722  {
723  return g;
724  }
725 
726  gidx_expand(gidx, fdistance);
727 
728  return gserialized_set_gidx(g, gidx);
729 }
730 
731 /***********************************************************************
732 * GiST N-D Index Operator Functions
733 */
734 
735 /*
736 * Do centroid to centroid n-d distance if you don't have
737 * re-check available (PgSQL 9.5+), do "real" n-d distance
738 * if you do
739 */
741 Datum gserialized_distance_nd(PG_FUNCTION_ARGS)
742 {
743  char b1mem[GIDX_MAX_SIZE];
744  GIDX *b1 = (GIDX*)b1mem;
745  char b2mem[GIDX_MAX_SIZE];
746  GIDX *b2 = (GIDX*)b2mem;
747 
748 #if POSTGIS_PGSQL_VERSION < 95
749 
750  /* Centroid-to-centroid distance */
751  Datum gs1 = PG_GETARG_DATUM(0);
752  Datum gs2 = PG_GETARG_DATUM(1);
753  double box_distance = FLT_MAX;
754 
755  /* Must be able to build box for each argument (ie, not empty geometry). */
756  if ( (gserialized_datum_get_gidx_p(gs1, b1) == LW_SUCCESS) &&
757  (gserialized_datum_get_gidx_p(gs2, b2) == LW_SUCCESS) )
758  {
759  box_distance = gidx_distance_leaf_centroid(b1, b2);
760  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(b1), gidx_to_string(b2));
761  }
762  PG_RETURN_FLOAT8(box_distance);
763 
764 #else /* POSTGIS_PGSQL_VERSION >= 96 */
765 
766  /* Feature-to-feature distance */
767  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
768  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
769  LWGEOM *lw1 = lwgeom_from_gserialized(geom1);
770  LWGEOM *lw2 = lwgeom_from_gserialized(geom2);
771  LWGEOM *closest;
772  double distance;
773 
774 
775  /* Find an exact shortest line w/ the dimensions we support */
776  if ( lwgeom_has_z(lw1) && lwgeom_has_z(lw2) )
777  {
778  closest = lwgeom_closest_line_3d(lw1, lw2);
779  distance = lwgeom_length(closest);
780  }
781  else
782  {
783  closest = lwgeom_closest_line(lw1, lw2);
784  distance = lwgeom_length_2d(closest);
785  }
786 
787  /* Un-sqrt the distance so we can add extra terms */
788  distance = distance*distance;
789 
790  /* Can only add the M term if both objects have M */
791  if ( lwgeom_has_m(lw1) && lwgeom_has_m(lw2) )
792  {
793  double m1, m2;
794  int usebox = false;
795 
796  if ( lwgeom_get_type(lw1) == POINTTYPE )
797  {
798  POINT4D p;
799  lwpoint_getPoint4d_p((LWPOINT*)lw1, &p);
800  m1 = p.m;
801  }
802  else if ( lwgeom_get_type(lw1) == LINETYPE )
803  {
804  LWPOINT *lwp1 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 0);
805  m1 = lwgeom_interpolate_point(lw1, lwp1);
806  lwpoint_free(lwp1);
807  }
808  else
809  {
810  usebox = true;
811  }
812 
813  if ( lwgeom_get_type(lw2) == POINTTYPE )
814  {
815  POINT4D p;
816  lwpoint_getPoint4d_p((LWPOINT*)lw2, &p);
817  m2 = p.m;
818  }
819  else if ( lwgeom_get_type(lw2) == LINETYPE )
820  {
821  LWPOINT *lwp2 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 1);
822  m2 = lwgeom_interpolate_point(lw2, lwp2);
823  lwpoint_free(lwp2);
824  }
825  else
826  {
827  usebox = true;
828  }
829 
830  if ( usebox )
831  {
832  double d;
833  gserialized_get_gidx_p(geom1, b1);
834  gserialized_get_gidx_p(geom2, b2);
835  d = gidx_distance_m(b1, b2);
836  distance += d*d;
837  }
838  else
839  {
840  distance += (m2-m1)*(m2-m1);
841  }
842  }
843 
844  lwgeom_free(closest);
845 
846  PG_FREE_IF_COPY(geom1, 0);
847  PG_FREE_IF_COPY(geom2, 1);
848  PG_RETURN_FLOAT8(sqrt(distance));
849 #endif /* POSTGIS_PGSQL_VERSION >= 96 */
850 }
851 
852 /*
853 ** '~' and operator function. Based on two serialized return true if
854 ** the first is contained by the second.
855 */
857 Datum gserialized_within(PG_FUNCTION_ARGS)
858 {
859  if ( gserialized_datum_predicate(PG_GETARG_DATUM(1), PG_GETARG_DATUM(0), gidx_contains) == LW_TRUE )
860  {
861  PG_RETURN_BOOL(TRUE);
862  }
863 
864  PG_RETURN_BOOL(FALSE);
865 }
866 
867 #if POSTGIS_PGSQL_VERSION > 94
868 /*
869 ** '~' and operator function. Based on a GIDX and a serialized return true if
870 ** the first is contained by the second.
871 */
872 PG_FUNCTION_INFO_V1(gserialized_gidx_geom_within);
873 Datum gserialized_gidx_geom_within(PG_FUNCTION_ARGS)
874 {
875  GIDX *gidx = (GIDX *)PG_GETARG_POINTER(0);
876 
877  if ( gserialized_datum_predicate_geom_gidx(PG_GETARG_DATUM(1), gidx, gidx_contains) == LW_TRUE )
878  PG_RETURN_BOOL(TRUE);
879 
880  PG_RETURN_BOOL(FALSE);
881 }
882 
883 /*
884 ** '~' and operator function. Based on two GIDX return true if
885 ** the first is contained by the second.
886 */
887 PG_FUNCTION_INFO_V1(gserialized_gidx_gidx_within);
888 Datum gserialized_gidx_gidx_within(PG_FUNCTION_ARGS)
889 {
890  if ( gidx_contains((GIDX *)PG_GETARG_POINTER(1), (GIDX *)PG_GETARG_POINTER(0)))
891  PG_RETURN_BOOL(TRUE);
892 
893  PG_RETURN_BOOL(FALSE);
894 }
895 #endif
896 
897 /*
898 ** '@' and operator function. Based on two serialized return true if
899 ** the first contains the second.
900 */
902 Datum gserialized_contains(PG_FUNCTION_ARGS)
903 {
904  if ( gserialized_datum_predicate(PG_GETARG_DATUM(0),PG_GETARG_DATUM(1), gidx_contains) == LW_TRUE )
905  {
906  PG_RETURN_BOOL(TRUE);
907  }
908 
909  PG_RETURN_BOOL(FALSE);
910 }
911 
912 #if POSTGIS_PGSQL_VERSION > 94
913 /*
914 ** '@' and operator function. Based on a GIDX and a serialized return true if
915 ** the first contains the second.
916 */
917 PG_FUNCTION_INFO_V1(gserialized_gidx_geom_contains);
918 Datum gserialized_gidx_geom_contains(PG_FUNCTION_ARGS)
919 {
920  GIDX *gidx = (GIDX *)PG_GETARG_POINTER(0);
921 
922  if ( gserialized_datum_predicate_gidx_geom(gidx, PG_GETARG_DATUM(1), gidx_contains) == LW_TRUE )
923  PG_RETURN_BOOL(TRUE);
924 
925  PG_RETURN_BOOL(FALSE);
926 }
927 
928 /*
929 ** '@' and operator function. Based on two GIDX return true if
930 ** the first contains the second.
931 */
932 PG_FUNCTION_INFO_V1(gserialized_gidx_gidx_contains);
933 Datum gserialized_gidx_gidx_contains(PG_FUNCTION_ARGS)
934 {
935  if ( gidx_contains((GIDX *)PG_GETARG_POINTER(0), (GIDX *)PG_GETARG_POINTER(1)))
936  PG_RETURN_BOOL(TRUE);
937 
938  PG_RETURN_BOOL(FALSE);
939 }
940 
941 PG_FUNCTION_INFO_V1(gserialized_gidx_geom_same);
942 Datum gserialized_gidx_geom_same(PG_FUNCTION_ARGS)
943 {
944  GIDX *gidx = (GIDX *)PG_GETARG_POINTER(0);
945 
946  if ( gserialized_datum_predicate_gidx_geom(gidx, PG_GETARG_DATUM(1), gidx_equals) == LW_TRUE )
947  PG_RETURN_BOOL(TRUE);
948 
949  PG_RETURN_BOOL(FALSE);
950 }
951 
952 PG_FUNCTION_INFO_V1(gserialized_gidx_gidx_same);
953 Datum gserialized_gidx_gidx_same(PG_FUNCTION_ARGS)
954 {
955  if ( gidx_equals((GIDX *)PG_GETARG_POINTER(0), (GIDX *)PG_GETARG_POINTER(1)) )
956  PG_RETURN_BOOL(TRUE);
957 
958  PG_RETURN_BOOL(FALSE);
959 }
960 #endif
961 
962 /*
963 ** '&&' operator function. Based on two serialized return true if
964 ** they overlap and false otherwise.
965 */
967 Datum gserialized_overlaps(PG_FUNCTION_ARGS)
968 {
969  if ( gserialized_datum_predicate(PG_GETARG_DATUM(0),PG_GETARG_DATUM(1), gidx_overlaps) == LW_TRUE )
970  {
971  PG_RETURN_BOOL(TRUE);
972  }
973 
974  PG_RETURN_BOOL(FALSE);
975 }
976 
977 #if POSTGIS_PGSQL_VERSION > 94
978 /*
979  * This is the cross-operator for the geographies
980  */
981 PG_FUNCTION_INFO_V1(gserialized_gidx_geog_overlaps);
982 Datum gserialized_gidx_geog_overlaps(PG_FUNCTION_ARGS)
983 {
984  GIDX *gidx = (GIDX *)PG_GETARG_POINTER(0);
985 
986  if ( gserialized_datum_predicate_gidx_geom(gidx, PG_GETARG_DATUM(1), gidx_overlaps) == LW_TRUE )
987  PG_RETURN_BOOL(TRUE);
988 
989  PG_RETURN_BOOL(FALSE);
990 }
991 
992 PG_FUNCTION_INFO_V1(gserialized_gidx_geom_overlaps);
993 Datum gserialized_gidx_geom_overlaps(PG_FUNCTION_ARGS)
994 {
995  GIDX *gidx = (GIDX *)PG_GETARG_POINTER(0);
996 
997  if ( gserialized_datum_predicate_gidx_geom(gidx, PG_GETARG_DATUM(1), gidx_overlaps) == LW_TRUE )
998  PG_RETURN_BOOL(TRUE);
999 
1000  PG_RETURN_BOOL(FALSE);
1001 }
1002 
1003 PG_FUNCTION_INFO_V1(gserialized_gidx_gidx_overlaps);
1004 Datum gserialized_gidx_gidx_overlaps(PG_FUNCTION_ARGS)
1005 {
1006  if ( gidx_overlaps((GIDX *)PG_GETARG_POINTER(0), (GIDX *)PG_GETARG_POINTER(1)) )
1007  PG_RETURN_BOOL(TRUE);
1008 
1009  PG_RETURN_BOOL(FALSE);
1010 }
1011 #endif
1012 
1013 /***********************************************************************
1014 * GiST Index Support Functions
1015 */
1016 
1017 /*
1018 ** GiST support function. Given a geography, return a "compressed"
1019 ** version. In this case, we convert the geography into a geocentric
1020 ** bounding box. If the geography already has the box embedded in it
1021 ** we pull that out and hand it back.
1022 */
1024 Datum gserialized_gist_compress(PG_FUNCTION_ARGS)
1025 {
1026  GISTENTRY *entry_in = (GISTENTRY*)PG_GETARG_POINTER(0);
1027  GISTENTRY *entry_out = NULL;
1028  char gidxmem[GIDX_MAX_SIZE];
1029  GIDX *bbox_out = (GIDX*)gidxmem;
1030  int result = LW_SUCCESS;
1031  int i;
1032 
1033  POSTGIS_DEBUG(4, "[GIST] 'compress' function called");
1034 
1035  /*
1036  ** Not a leaf key? There's nothing to do.
1037  ** Return the input unchanged.
1038  */
1039  if ( ! entry_in->leafkey )
1040  {
1041  POSTGIS_DEBUG(4, "[GIST] non-leafkey entry, returning input unaltered");
1042  PG_RETURN_POINTER(entry_in);
1043  }
1044 
1045  POSTGIS_DEBUG(4, "[GIST] processing leafkey input");
1046  entry_out = palloc(sizeof(GISTENTRY));
1047 
1048  /*
1049  ** Null key? Make a copy of the input entry and
1050  ** return.
1051  */
1052  if ( DatumGetPointer(entry_in->key) == NULL )
1053  {
1054  POSTGIS_DEBUG(4, "[GIST] leafkey is null");
1055  gistentryinit(*entry_out, (Datum) 0, entry_in->rel,
1056  entry_in->page, entry_in->offset, FALSE);
1057  POSTGIS_DEBUG(4, "[GIST] returning copy of input");
1058  PG_RETURN_POINTER(entry_out);
1059  }
1060 
1061  /* Extract our index key from the GiST entry. */
1062  result = gserialized_datum_get_gidx_p(entry_in->key, bbox_out);
1063 
1064  /* Is the bounding box valid (non-empty, non-infinite) ?
1065  * If not, use the "unknown" GIDX. */
1066  if ( result == LW_FAILURE )
1067  {
1068  POSTGIS_DEBUG(4, "[GIST] empty geometry!");
1069  gidx_set_unknown(bbox_out);
1070  gistentryinit(*entry_out, PointerGetDatum(gidx_copy(bbox_out)),
1071  entry_in->rel, entry_in->page,
1072  entry_in->offset, FALSE);
1073  PG_RETURN_POINTER(entry_out);
1074  }
1075 
1076  POSTGIS_DEBUGF(4, "[GIST] got entry_in->key: %s", gidx_to_string(bbox_out));
1077 
1078  /* Check all the dimensions for finite values.
1079  * If not, use the "unknown" GIDX as a key */
1080  for ( i = 0; i < GIDX_NDIMS(bbox_out); i++ )
1081  {
1082  if ( ! isfinite(GIDX_GET_MAX(bbox_out, i))
1083  || ! isfinite(GIDX_GET_MIN(bbox_out, i)) )
1084  {
1085  gidx_set_unknown(bbox_out);
1086  gistentryinit(*entry_out,
1087  PointerGetDatum(gidx_copy(bbox_out)),
1088  entry_in->rel, entry_in->page,
1089  entry_in->offset, FALSE);
1090  PG_RETURN_POINTER(entry_out);
1091  }
1092  }
1093 
1094  /* Enure bounding box has minimums below maximums. */
1095  gidx_validate(bbox_out);
1096 
1097  /* Prepare GISTENTRY for return. */
1098  gistentryinit(*entry_out, PointerGetDatum(gidx_copy(bbox_out)),
1099  entry_in->rel, entry_in->page, entry_in->offset, FALSE);
1100 
1101  /* Return GISTENTRY. */
1102  POSTGIS_DEBUG(4, "[GIST] 'compress' function complete");
1103  PG_RETURN_POINTER(entry_out);
1104 }
1105 
1106 /*
1107 ** GiST support function.
1108 ** Decompress an entry. Unused for geography, so we return.
1109 */
1111 Datum gserialized_gist_decompress(PG_FUNCTION_ARGS)
1112 {
1113  POSTGIS_DEBUG(5, "[GIST] 'decompress' function called");
1114  /* We don't decompress. Just return the input. */
1115  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
1116 }
1117 
1118 /*
1119 ** GiST support function. Called from gserialized_gist_consistent below.
1120 */
1121 static inline bool gserialized_gist_consistent_leaf(GIDX *key, GIDX *query, StrategyNumber strategy)
1122 {
1123  bool retval;
1124 
1125  POSTGIS_DEBUGF(4, "[GIST] leaf consistent, strategy [%d], count[%i]",
1126  strategy, geog_counter_leaf++);
1127 
1128  switch (strategy)
1129  {
1130  case RTOverlapStrategyNumber:
1131  retval = (bool) gidx_overlaps(key, query);
1132  break;
1133  case RTSameStrategyNumber:
1134  retval = (bool) gidx_equals(key, query);
1135  break;
1136  case RTContainsStrategyNumber:
1137  case RTOldContainsStrategyNumber:
1138  retval = (bool) gidx_contains(key, query);
1139  break;
1140  case RTContainedByStrategyNumber:
1141  case RTOldContainedByStrategyNumber:
1142  retval = (bool) gidx_contains(query, key);
1143  break;
1144  default:
1145  retval = FALSE;
1146  }
1147 
1148  return (retval);
1149 }
1150 
1151 /*
1152 ** GiST support function. Called from gserialized_gist_consistent below.
1153 */
1154 static inline bool gserialized_gist_consistent_internal(GIDX *key, GIDX *query, StrategyNumber strategy)
1155 {
1156  bool retval;
1157 
1158  POSTGIS_DEBUGF(4, "[GIST] internal consistent, strategy [%d], count[%i], query[%s], key[%s]",
1159  strategy, geog_counter_internal++, gidx_to_string(query), gidx_to_string(key) );
1160 
1161  switch (strategy)
1162  {
1163  case RTOverlapStrategyNumber:
1164  retval = (bool) gidx_overlaps(key, query);
1165  break;
1166  case RTSameStrategyNumber:
1167  case RTContainsStrategyNumber:
1168  case RTOldContainsStrategyNumber:
1169  retval = (bool) gidx_contains(key, query);
1170  break;
1171  case RTContainedByStrategyNumber:
1172  case RTOldContainedByStrategyNumber:
1173  retval = (bool) gidx_overlaps(key, query);
1174  break;
1175  default:
1176  retval = FALSE;
1177  }
1178 
1179  return (retval);
1180 }
1181 
1182 /*
1183 ** GiST support function. Take in a query and an entry and see what the
1184 ** relationship is, based on the query strategy.
1185 */
1187 Datum gserialized_gist_consistent(PG_FUNCTION_ARGS)
1188 {
1189  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1190  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1191  bool result;
1192  char gidxmem[GIDX_MAX_SIZE];
1193  GIDX *query_gbox_index = (GIDX*)gidxmem;
1194 
1195  /* PostgreSQL 8.4 and later require the RECHECK flag to be set here,
1196  rather than being supplied as part of the operator class definition */
1197  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1198 
1199  /* We set recheck to false to avoid repeatedly pulling every "possibly matched" geometry
1200  out during index scans. For cases when the geometries are large, rechecking
1201  can make things twice as slow. */
1202  *recheck = false;
1203 
1204  POSTGIS_DEBUG(4, "[GIST] 'consistent' function called");
1205 
1206  /* Quick sanity check on query argument. */
1207  if ( DatumGetPointer(PG_GETARG_DATUM(1)) == NULL )
1208  {
1209  POSTGIS_DEBUG(4, "[GIST] null query pointer (!?!), returning false");
1210  PG_RETURN_BOOL(FALSE); /* NULL query! This is screwy! */
1211  }
1212 
1213  /* Quick sanity check on entry key. */
1214  if ( DatumGetPointer(entry->key) == NULL )
1215  {
1216  POSTGIS_DEBUG(4, "[GIST] null index entry, returning false");
1217  PG_RETURN_BOOL(FALSE); /* NULL entry! */
1218  }
1219 
1220  /* Null box should never make this far. */
1221  if ( gserialized_datum_get_gidx_p(PG_GETARG_DATUM(1), query_gbox_index) == LW_FAILURE )
1222  {
1223  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1224  PG_RETURN_BOOL(FALSE);
1225  }
1226 
1227  /* Treat leaf node tests different from internal nodes */
1228  if (GIST_LEAF(entry))
1229  {
1231  (GIDX*)DatumGetPointer(entry->key),
1232  query_gbox_index, strategy);
1233  }
1234  else
1235  {
1237  (GIDX*)DatumGetPointer(entry->key),
1238  query_gbox_index, strategy);
1239  }
1240 
1241  PG_RETURN_BOOL(result);
1242 }
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(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  GIDX *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 = (GIDX*)DatumGetPointer(origentry->key);
1262  gbox_index_new = (GIDX*)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 (volume difference in this case). */
1273  size_union = gidx_union_volume(gbox_index_orig, gbox_index_new);
1274  size_orig = gidx_volume(gbox_index_orig);
1275  *result = size_union - size_orig;
1276 
1277  /* All things being equal, we prefer to expand small boxes rather than large boxes.
1278  NOTE: This code seemed to be causing badly balanced trees to be built
1279  and therefore has been commented out. Not sure why it didn't work,
1280  but it didn't.
1281  if( FP_IS_ZERO(*result) )
1282  if( FP_IS_ZERO(size_orig) )
1283  *result = 0.0;
1284  else
1285  *result = 1.0 - (1.0/(1.0 + size_orig));
1286  else
1287  *result += 1.0;
1288  */
1289 
1290  POSTGIS_DEBUGF(4, "[GIST] union size (%.12f), original size (%.12f), penalty (%.12f)", size_union, size_orig, *result);
1291 
1292  PG_RETURN_POINTER(result);
1293 }
1294 
1295 /*
1296 ** GiST support function. Merge all the boxes in a page into one master box.
1297 */
1299 Datum gserialized_gist_union(PG_FUNCTION_ARGS)
1300 {
1301  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
1302  int *sizep = (int *) PG_GETARG_POINTER(1); /* Size of the return value */
1303  int numranges, i;
1304  GIDX *box_cur, *box_union;
1305 
1306  POSTGIS_DEBUG(4, "[GIST] 'union' function called");
1307 
1308  numranges = entryvec->n;
1309 
1310  box_cur = (GIDX*) DatumGetPointer(entryvec->vector[0].key);
1311 
1312  box_union = gidx_copy(box_cur);
1313 
1314  for ( i = 1; i < numranges; i++ )
1315  {
1316  box_cur = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
1317  gidx_merge(&box_union, box_cur);
1318  }
1319 
1320  *sizep = VARSIZE(box_union);
1321 
1322  POSTGIS_DEBUGF(4, "[GIST] union called, numranges(%i), pageunion %s", numranges, gidx_to_string(box_union));
1323 
1324  PG_RETURN_POINTER(box_union);
1325 
1326 }
1327 
1328 /*
1329 ** GiST support function. Test equality of keys.
1330 */
1332 Datum gserialized_gist_same(PG_FUNCTION_ARGS)
1333 {
1334  GIDX *b1 = (GIDX*)PG_GETARG_POINTER(0);
1335  GIDX *b2 = (GIDX*)PG_GETARG_POINTER(1);
1336  bool *result = (bool*)PG_GETARG_POINTER(2);
1337 
1338  POSTGIS_DEBUG(4, "[GIST] 'same' function called");
1339 
1340  *result = gidx_equals(b1, b2);
1341 
1342  PG_RETURN_POINTER(result);
1343 }
1344 
1345 
1346 
1347 
1349 Datum gserialized_gist_geog_distance(PG_FUNCTION_ARGS)
1350 {
1351  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1352  Datum query_datum = PG_GETARG_DATUM(1);
1353  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1354 #if POSTGIS_PGSQL_VERSION >= 95
1355  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1356 #endif
1357  char query_box_mem[GIDX_MAX_SIZE];
1358  GIDX *query_box = (GIDX*)query_box_mem;
1359  GIDX *entry_box;
1360  double distance;
1361 
1362  POSTGIS_DEBUGF(3, "[GIST] '%s' function called", __func__);
1363 
1364  /* We are using '13' as the gist geography distance <-> strategy number */
1365  if ( strategy != 13 )
1366  {
1367  elog(ERROR, "unrecognized strategy number: %d", strategy);
1368  PG_RETURN_FLOAT8(FLT_MAX);
1369  }
1370 
1371  /* Null box should never make this far. */
1372  if ( gserialized_datum_get_gidx_p(query_datum, query_box) == LW_FAILURE )
1373  {
1374  POSTGIS_DEBUG(2, "[GIST] null query_gbox_index!");
1375  PG_RETURN_FLOAT8(FLT_MAX);
1376  }
1377 
1378 #if POSTGIS_PGSQL_VERSION >= 95
1379  /* When we hit leaf nodes, it's time to turn on recheck */
1380  if (GIST_LEAF(entry))
1381  {
1382  *recheck = true;
1383  }
1384 #endif
1385 
1386  /* Get the entry box */
1387  entry_box = (GIDX*)DatumGetPointer(entry->key);
1388 
1389  /* Return distances from key-based tests should always be */
1390  /* the minimum possible distance, box-to-box */
1391  /* We scale up to "world units" so that the box-to-box distances */
1392  /* compare reasonably with the over-the-spheroid distances that */
1393  /* the recheck process will turn up */
1394  distance = WGS84_RADIUS * gidx_distance(entry_box, query_box, 0);
1395  POSTGIS_DEBUGF(2, "[GIST] '%s' got distance %g", __func__, distance);
1396 
1397  PG_RETURN_FLOAT8(distance);
1398 }
1399 
1400 
1401 /*
1402 ** GiST support function.
1403 ** Take in a query and an entry and return the "distance" between them.
1404 **
1405 ** Given an index entry p and a query value q, this function determines the
1406 ** index entry's "distance" from the query value. This function must be
1407 ** supplied if the operator class contains any ordering operators. A query
1408 ** using the ordering operator will be implemented by returning index entries
1409 ** with the smallest "distance" values first, so the results must be consistent
1410 ** with the operator's semantics. For a leaf index entry the result just
1411 ** represents the distance to the index entry; for an internal tree node, the
1412 ** result must be the smallest distance that any child entry could have.
1413 **
1414 ** Strategy 13 is centroid-based distance tests <<->>
1415 ** Strategy 20 is cpa based distance tests |=|
1416 */
1418 Datum gserialized_gist_distance(PG_FUNCTION_ARGS)
1419 {
1420  GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
1421  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1422  char query_box_mem[GIDX_MAX_SIZE];
1423  GIDX *query_box = (GIDX*)query_box_mem;
1424  GIDX *entry_box;
1425 #if POSTGIS_PGSQL_VERSION >= 95
1426  bool *recheck = (bool *) PG_GETARG_POINTER(4);
1427 #endif
1428 
1429  double distance;
1430 
1431  POSTGIS_DEBUG(4, "[GIST] 'distance' function called");
1432 
1433  /* Strategy 13 is <<->> */
1434  /* Strategy 20 is |=| */
1435  if ( strategy != 13 && strategy != 20 ) {
1436  elog(ERROR, "unrecognized strategy number: %d", strategy);
1437  PG_RETURN_FLOAT8(FLT_MAX);
1438  }
1439 
1440  /* Null box should never make this far. */
1441  if ( gserialized_datum_get_gidx_p(PG_GETARG_DATUM(1), query_box) == LW_FAILURE )
1442  {
1443  POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
1444  PG_RETURN_FLOAT8(FLT_MAX);
1445  }
1446 
1447  /* Get the entry box */
1448  entry_box = (GIDX*)DatumGetPointer(entry->key);
1449 
1450 #if POSTGIS_PGSQL_VERSION >= 95
1451 
1452  /* Strategy 20 is |=| */
1453  distance = gidx_distance(entry_box, query_box, strategy == 20);
1454 
1455  /* Treat leaf node tests different from internal nodes */
1456  if (GIST_LEAF(entry))
1457  *recheck = true;
1458 #else
1459 
1460  if ( strategy == 20 )
1461  {
1462  elog(ERROR, "You need PostgreSQL 9.5.0 or higher in order to use |=| with index");
1463  PG_RETURN_FLOAT8(FLT_MAX);
1464  }
1465 
1466  /* Treat leaf node tests different from internal nodes */
1467  if (GIST_LEAF(entry))
1468  {
1469  /* Calculate distance to leaves */
1470  distance = (double)gidx_distance_leaf_centroid(entry_box, query_box);
1471  }
1472  else
1473  {
1474  /* Calculate distance for internal nodes */
1475  distance = (double)gidx_distance_node_centroid(entry_box, query_box);
1476  }
1477 #endif
1478  PG_RETURN_FLOAT8(distance);
1479 }
1480 
1481 
1482 /*
1483 ** Utility function to add entries to the axis partition lists in the
1484 ** picksplit function.
1485 */
1486 static void gserialized_gist_picksplit_addlist(OffsetNumber *list, GIDX **box_union, GIDX *box_current, int *pos, int num)
1487 {
1488  if ( *pos )
1489  gidx_merge(box_union, box_current);
1490  else
1491  memcpy((void*)(*box_union), (void*)box_current, VARSIZE(box_current));
1492  list[*pos] = num;
1493  (*pos)++;
1494 }
1495 
1496 /*
1497 ** Utility function check whether the number of entries two halves of the
1498 ** space constitute a "bad ratio" (poor balance).
1499 */
1501 {
1502  POSTGIS_DEBUGF(4, "[GIST] checking split ratio (%d, %d)", x, y);
1503  if ( (y == 0) || (((float)x / (float)y) < LIMIT_RATIO) ||
1504  (x == 0) || (((float)y / (float)x) < LIMIT_RATIO) )
1505  return TRUE;
1506 
1507  return FALSE;
1508 }
1509 
1510 static bool gserialized_gist_picksplit_badratios(int *pos, int dims)
1511 {
1512  int i;
1513  for ( i = 0; i < dims; i++ )
1514  {
1515  if ( gserialized_gist_picksplit_badratio(pos[2*i],pos[2*i+1]) == FALSE )
1516  return FALSE;
1517  }
1518  return TRUE;
1519 }
1520 
1521 
1522 /*
1523 ** Where the picksplit algorithm cannot find any basis for splitting one way
1524 ** or another, we simply split the overflowing node in half.
1525 */
1526 static void gserialized_gist_picksplit_fallback(GistEntryVector *entryvec, GIST_SPLITVEC *v)
1527 {
1528  OffsetNumber i, maxoff;
1529  GIDX *unionL = NULL;
1530  GIDX *unionR = NULL;
1531  int nbytes;
1532 
1533  POSTGIS_DEBUG(4, "[GIST] in fallback picksplit function");
1534 
1535  maxoff = entryvec->n - 1;
1536 
1537  nbytes = (maxoff + 2) * sizeof(OffsetNumber);
1538  v->spl_left = (OffsetNumber*) palloc(nbytes);
1539  v->spl_right = (OffsetNumber*) palloc(nbytes);
1540  v->spl_nleft = v->spl_nright = 0;
1541 
1542  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
1543  {
1544  GIDX *cur = (GIDX*)DatumGetPointer(entryvec->vector[i].key);
1545 
1546  if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
1547  {
1548  v->spl_left[v->spl_nleft] = i;
1549  if (unionL == NULL)
1550  {
1551  unionL = gidx_copy(cur);
1552  }
1553  else
1554  {
1555  gidx_merge(&unionL, cur);
1556  }
1557  v->spl_nleft++;
1558  }
1559  else
1560  {
1561  v->spl_right[v->spl_nright] = i;
1562  if (unionR == NULL)
1563  {
1564  unionR = gidx_copy(cur);
1565  }
1566  else
1567  {
1568  gidx_merge(&unionR, cur);
1569  }
1570  v->spl_nright++;
1571  }
1572  }
1573 
1574  if (v->spl_ldatum_exists)
1575  gidx_merge(&unionL, (GIDX*)DatumGetPointer(v->spl_ldatum));
1576 
1577  v->spl_ldatum = PointerGetDatum(unionL);
1578 
1579  if (v->spl_rdatum_exists)
1580  gidx_merge(&unionR, (GIDX*)DatumGetPointer(v->spl_rdatum));
1581 
1582  v->spl_rdatum = PointerGetDatum(unionR);
1583  v->spl_ldatum_exists = v->spl_rdatum_exists = false;
1584 }
1585 
1586 
1587 
1588 static void gserialized_gist_picksplit_constructsplit(GIST_SPLITVEC *v, OffsetNumber *list1, int nlist1, GIDX **union1, OffsetNumber *list2, int nlist2, GIDX **union2)
1589 {
1590  bool firstToLeft = true;
1591 
1592  POSTGIS_DEBUG(4, "[GIST] picksplit in constructsplit function");
1593 
1594  if (v->spl_ldatum_exists || v->spl_rdatum_exists)
1595  {
1596  if (v->spl_ldatum_exists && v->spl_rdatum_exists)
1597  {
1598  GIDX *LRl = gidx_copy(*union1);
1599  GIDX *LRr = gidx_copy(*union2);
1600  GIDX *RLl = gidx_copy(*union2);
1601  GIDX *RLr = gidx_copy(*union1);
1602  double sizeLR, sizeRL;
1603 
1604  gidx_merge(&LRl, (GIDX*)DatumGetPointer(v->spl_ldatum));
1605  gidx_merge(&LRr, (GIDX*)DatumGetPointer(v->spl_rdatum));
1606  gidx_merge(&RLl, (GIDX*)DatumGetPointer(v->spl_ldatum));
1607  gidx_merge(&RLr, (GIDX*)DatumGetPointer(v->spl_rdatum));
1608 
1609  sizeLR = gidx_inter_volume(LRl,LRr);
1610  sizeRL = gidx_inter_volume(RLl,RLr);
1611 
1612  POSTGIS_DEBUGF(4, "[GIST] sizeLR / sizeRL == %.12g / %.12g", sizeLR, sizeRL);
1613 
1614  if (sizeLR > sizeRL)
1615  firstToLeft = false;
1616 
1617  }
1618  else
1619  {
1620  float p1, p2;
1621  GISTENTRY oldUnion, addon;
1622 
1623  gistentryinit(oldUnion, (v->spl_ldatum_exists) ? v->spl_ldatum : v->spl_rdatum,
1624  NULL, NULL, InvalidOffsetNumber, FALSE);
1625 
1626  gistentryinit(addon, PointerGetDatum(*union1), NULL, NULL, InvalidOffsetNumber, FALSE);
1627  DirectFunctionCall3(gserialized_gist_penalty, PointerGetDatum(&oldUnion), PointerGetDatum(&addon), PointerGetDatum(&p1));
1628  gistentryinit(addon, PointerGetDatum(*union2), NULL, NULL, InvalidOffsetNumber, FALSE);
1629  DirectFunctionCall3(gserialized_gist_penalty, PointerGetDatum(&oldUnion), PointerGetDatum(&addon), PointerGetDatum(&p2));
1630 
1631  POSTGIS_DEBUGF(4, "[GIST] p1 / p2 == %.12g / %.12g", p1, p2);
1632 
1633  if ((v->spl_ldatum_exists && p1 > p2) || (v->spl_rdatum_exists && p1 < p2))
1634  firstToLeft = false;
1635  }
1636  }
1637 
1638  POSTGIS_DEBUGF(4, "[GIST] firstToLeft == %d", firstToLeft);
1639 
1640  if (firstToLeft)
1641  {
1642  v->spl_left = list1;
1643  v->spl_right = list2;
1644  v->spl_nleft = nlist1;
1645  v->spl_nright = nlist2;
1646  if (v->spl_ldatum_exists)
1647  gidx_merge(union1, (GIDX*)DatumGetPointer(v->spl_ldatum));
1648  v->spl_ldatum = PointerGetDatum(*union1);
1649  if (v->spl_rdatum_exists)
1650  gidx_merge(union2, (GIDX*)DatumGetPointer(v->spl_rdatum));
1651  v->spl_rdatum = PointerGetDatum(*union2);
1652  }
1653  else
1654  {
1655  v->spl_left = list2;
1656  v->spl_right = list1;
1657  v->spl_nleft = nlist2;
1658  v->spl_nright = nlist1;
1659  if (v->spl_ldatum_exists)
1660  gidx_merge(union2, (GIDX*)DatumGetPointer(v->spl_ldatum));
1661  v->spl_ldatum = PointerGetDatum(*union2);
1662  if (v->spl_rdatum_exists)
1663  gidx_merge(union1, (GIDX*)DatumGetPointer(v->spl_rdatum));
1664  v->spl_rdatum = PointerGetDatum(*union1);
1665  }
1666 
1667  v->spl_ldatum_exists = v->spl_rdatum_exists = false;
1668 }
1669 
1670 
1671 #define BELOW(d) (2*(d))
1672 #define ABOVE(d) ((2*(d))+1)
1673 
1674 /*
1675 ** GiST support function. Split an overflowing node into two new nodes.
1676 ** Uses linear algorithm from Ang & Tan [2], dividing node extent into
1677 ** four quadrants, and splitting along the axis that most evenly distributes
1678 ** entries between the new nodes.
1679 ** TODO: Re-evaluate this in light of R*Tree picksplit approaches.
1680 */
1682 Datum gserialized_gist_picksplit(PG_FUNCTION_ARGS)
1683 {
1684 
1685  GistEntryVector *entryvec = (GistEntryVector*) PG_GETARG_POINTER(0);
1686 
1687  GIST_SPLITVEC *v = (GIST_SPLITVEC*) PG_GETARG_POINTER(1);
1688  OffsetNumber i;
1689  /* One union box for each half of the space. */
1690  GIDX **box_union;
1691  /* One offset number list for each half of the space. */
1692  OffsetNumber **list;
1693  /* One position index for each half of the space. */
1694  int *pos;
1695  GIDX *box_pageunion;
1696  GIDX *box_current;
1697  int direction = -1;
1698  bool all_entries_equal = true;
1699  OffsetNumber max_offset;
1700  int nbytes, ndims_pageunion, d;
1701  int posmin = entryvec->n;
1702 
1703  POSTGIS_DEBUG(4, "[GIST] 'picksplit' function called");
1704 
1705  /*
1706  ** First calculate the bounding box and maximum number of dimensions in this page.
1707  */
1708 
1709  max_offset = entryvec->n - 1;
1710  box_current = (GIDX*) DatumGetPointer(entryvec->vector[FirstOffsetNumber].key);
1711  box_pageunion = gidx_copy(box_current);
1712 
1713  /* Calculate the containing box (box_pageunion) for the whole page we are going to split. */
1714  for ( i = OffsetNumberNext(FirstOffsetNumber); i <= max_offset; i = OffsetNumberNext(i) )
1715  {
1716  box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
1717 
1718  if ( all_entries_equal == true && ! gidx_equals (box_pageunion, box_current) )
1719  all_entries_equal = false;
1720 
1721  gidx_merge( &box_pageunion, box_current );
1722  }
1723 
1724  POSTGIS_DEBUGF(3, "[GIST] box_pageunion: %s", gidx_to_string(box_pageunion));
1725 
1726  /* Every box in the page is the same! So, we split and just put half the boxes in each child. */
1727  if ( all_entries_equal )
1728  {
1729  POSTGIS_DEBUG(4, "[GIST] picksplit finds all entries equal!");
1731  PG_RETURN_POINTER(v);
1732  }
1733 
1734  /* Initialize memory structures. */
1735  nbytes = (max_offset + 2) * sizeof(OffsetNumber);
1736  ndims_pageunion = GIDX_NDIMS(box_pageunion);
1737  POSTGIS_DEBUGF(4, "[GIST] ndims_pageunion == %d", ndims_pageunion);
1738  pos = palloc(2*ndims_pageunion * sizeof(int));
1739  list = palloc(2*ndims_pageunion * sizeof(OffsetNumber*));
1740  box_union = palloc(2*ndims_pageunion * sizeof(GIDX*));
1741  for ( d = 0; d < ndims_pageunion; d++ )
1742  {
1743  list[BELOW(d)] = (OffsetNumber*) palloc(nbytes);
1744  list[ABOVE(d)] = (OffsetNumber*) palloc(nbytes);
1745  box_union[BELOW(d)] = gidx_new(ndims_pageunion);
1746  box_union[ABOVE(d)] = gidx_new(ndims_pageunion);
1747  pos[BELOW(d)] = 0;
1748  pos[ABOVE(d)] = 0;
1749  }
1750 
1751  /*
1752  ** Assign each entry in the node to the volume partitions it belongs to,
1753  ** such as "above the x/y plane, left of the y/z plane, below the x/z plane".
1754  ** Each entry thereby ends up in three of the six partitions.
1755  */
1756  POSTGIS_DEBUG(4, "[GIST] 'picksplit' calculating best split axis");
1757  for ( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
1758  {
1759  box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
1760 
1761  for ( d = 0; d < ndims_pageunion; d++ )
1762  {
1763  if ( GIDX_GET_MIN(box_current,d)-GIDX_GET_MIN(box_pageunion,d) < GIDX_GET_MAX(box_pageunion,d)-GIDX_GET_MAX(box_current,d) )
1764  {
1765  gserialized_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
1766  }
1767  else
1768  {
1769  gserialized_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
1770  }
1771 
1772  }
1773 
1774  }
1775 
1776  /*
1777  ** "Bad disposition", too many entries fell into one octant of the space, so no matter which
1778  ** plane we choose to split on, we're going to end up with a mostly full node. Where the
1779  ** data is pretty homogeneous (lots of duplicates) entries that are equidistant from the
1780  ** sides of the page union box can occasionally all end up in one place, leading
1781  ** to this condition.
1782  */
1783  if ( gserialized_gist_picksplit_badratios(pos,ndims_pageunion) == TRUE )
1784  {
1785  /*
1786  ** Instead we split on center points and see if we do better.
1787  ** First calculate the average center point for each axis.
1788  */
1789  double *avgCenter = palloc(ndims_pageunion * sizeof(double));
1790 
1791  for ( d = 0; d < ndims_pageunion; d++ )
1792  {
1793  avgCenter[d] = 0.0;
1794  }
1795 
1796  POSTGIS_DEBUG(4, "[GIST] picksplit can't find good split axis, trying center point method");
1797 
1798  for ( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
1799  {
1800  box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
1801  for ( d = 0; d < ndims_pageunion; d++ )
1802  {
1803  avgCenter[d] += (GIDX_GET_MAX(box_current,d) + GIDX_GET_MIN(box_current,d)) / 2.0;
1804  }
1805  }
1806  for ( d = 0; d < ndims_pageunion; d++ )
1807  {
1808  avgCenter[d] /= max_offset;
1809  pos[BELOW(d)] = pos[ABOVE(d)] = 0; /* Re-initialize our counters. */
1810  POSTGIS_DEBUGF(4, "[GIST] picksplit average center point[%d] = %.12g", d, avgCenter[d]);
1811  }
1812 
1813  /* For each of our entries... */
1814  for ( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
1815  {
1816  double center;
1817  box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
1818 
1819  for ( d = 0; d < ndims_pageunion; d++ )
1820  {
1821  center = (GIDX_GET_MIN(box_current,d)+GIDX_GET_MAX(box_current,d))/2.0;
1822  if ( center < avgCenter[d] )
1823  gserialized_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
1824  else if ( FPeq(center, avgCenter[d]) )
1825  if ( pos[BELOW(d)] > pos[ABOVE(d)] )
1826  gserialized_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
1827  else
1828  gserialized_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
1829  else
1830  gserialized_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
1831  }
1832 
1833  }
1834 
1835  /* Do we have a good disposition now? If not, screw it, just cut the node in half. */
1836  if ( gserialized_gist_picksplit_badratios(pos,ndims_pageunion) == TRUE )
1837  {
1838  POSTGIS_DEBUG(4, "[GIST] picksplit still cannot find a good split! just cutting the node in half");
1840  PG_RETURN_POINTER(v);
1841  }
1842 
1843  }
1844 
1845  /*
1846  ** Now, what splitting plane gives us the most even ratio of
1847  ** entries in our child pages? Since each split region has been apportioned entries
1848  ** against the same number of total entries, the axis that has the smallest maximum
1849  ** number of entries in its regions is the most evenly distributed.
1850  ** TODO: what if the distributions are equal in two or more axes?
1851  */
1852  for ( d = 0; d < ndims_pageunion; d++ )
1853  {
1854  int posd = Max(pos[ABOVE(d)],pos[BELOW(d)]);
1855  if ( posd < posmin )
1856  {
1857  direction = d;
1858  posmin = posd;
1859  }
1860  }
1861  if ( direction == -1 || posmin == entryvec->n )
1862  {
1863  /* ERROR OUT HERE */
1864  elog(ERROR, "Error in building split, unable to determine split direction.");
1865  }
1866 
1867  POSTGIS_DEBUGF(3, "[GIST] 'picksplit' splitting on axis %d", direction);
1868 
1870  pos[BELOW(direction)],
1871  &(box_union[BELOW(direction)]),
1872  list[ABOVE(direction)],
1873  pos[ABOVE(direction)],
1874  &(box_union[ABOVE(direction)]) );
1875 
1876  POSTGIS_DEBUGF(4, "[GIST] spl_ldatum: %s", gidx_to_string((GIDX*)v->spl_ldatum));
1877  POSTGIS_DEBUGF(4, "[GIST] spl_rdatum: %s", gidx_to_string((GIDX*)v->spl_rdatum));
1878 
1879  POSTGIS_DEBUGF(4, "[GIST] axis %d: parent range (%.12g, %.12g) left range (%.12g, %.12g), right range (%.12g, %.12g)",
1880  direction,
1881  GIDX_GET_MIN(box_pageunion, direction), GIDX_GET_MAX(box_pageunion, direction),
1882  GIDX_GET_MIN((GIDX*)v->spl_ldatum, direction), GIDX_GET_MAX((GIDX*)v->spl_ldatum, direction),
1883  GIDX_GET_MIN((GIDX*)v->spl_rdatum, direction), GIDX_GET_MAX((GIDX*)v->spl_rdatum, direction) );
1884 
1885  PG_RETURN_POINTER(v);
1886 
1887 }
1888 
1889 /*
1890 ** The GIDX key must be defined as a PostgreSQL type, even though it is only
1891 ** ever used internally. These no-op stubs are used to bind the type.
1892 */
1894 Datum gidx_in(PG_FUNCTION_ARGS)
1895 {
1896  ereport(ERROR,(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1897  errmsg("function gidx_in not implemented")));
1898  PG_RETURN_POINTER(NULL);
1899 }
1900 
1902 Datum gidx_out(PG_FUNCTION_ARGS)
1903 {
1904  GIDX *box = (GIDX *) PG_GETARG_POINTER(0);
1905  char *result = gidx_to_string(box);
1906  PG_RETURN_CSTRING(result);
1907 }
#define LINETYPE
Definition: liblwgeom.h:85
static void gidx_merge(GIDX **b_union, GIDX *b_new)
static GIDX * gidx_copy(GIDX *b)
static double gidx_distance_node_centroid(const GIDX *node, const GIDX *query)
Datum gserialized_within(PG_FUNCTION_ARGS)
double m
Definition: liblwgeom.h:351
Datum gserialized_gist_decompress(PG_FUNCTION_ARGS)
static bool gidx_is_unknown(const GIDX *a)
Datum gserialized_gist_geog_distance(PG_FUNCTION_ARGS)
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
#define WGS84_RADIUS
Definition: liblwgeom.h:130
PG_FUNCTION_INFO_V1(gserialized_distance_nd)
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:842
static float gidx_volume(GIDX *a)
Datum gserialized_contains(PG_FUNCTION_ARGS)
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:195
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
static void gserialized_gist_picksplit_constructsplit(GIST_SPLITVEC *v, OffsetNumber *list1, int nlist1, GIDX **union1, OffsetNumber *list2, int nlist2, GIDX **union2)
Datum gserialized_gist_penalty(PG_FUNCTION_ARGS)
#define LW_SUCCESS
Definition: liblwgeom.h:79
static void gidx_set_unknown(GIDX *a)
static float gidx_union_volume(GIDX *a, GIDX *b)
LWGEOM * lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:88
bool gidx_contains(GIDX *a, GIDX *b)
static int gserialized_datum_predicate(Datum gs1, Datum gs2, gidx_predicate predicate)
Support function.
double lwgeom_length_2d(const LWGEOM *geom)
Definition: lwgeom.c:1681
Datum gserialized_gist_union(PG_FUNCTION_ARGS)
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
static bool gserialized_gist_consistent_leaf(GIDX *key, GIDX *query, StrategyNumber strategy)
double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
Find the measure value at the location on the line closest to the point.
#define LW_FAILURE
Definition: liblwgeom.h:78
static void gidx_dimensionality_check(GIDX **a, GIDX **b)
static bool gidx_equals(GIDX *a, GIDX *b)
GSERIALIZED * gserialized_expand(GSERIALIZED *g, double distance)
Return a GSERIALIZED with an expanded bounding box.
#define LIMIT_RATIO
static void gidx_validate(GIDX *b)
Datum gserialized_gist_distance(PG_FUNCTION_ARGS)
Datum gserialized_gist_compress(PG_FUNCTION_ARGS)
Datum gserialized_distance_nd(PG_FUNCTION_ARGS)
#define LW_FALSE
Definition: liblwgeom.h:76
#define BELOW(d)
static double gidx_distance(const GIDX *a, const GIDX *b, int m_is_time)
Calculate the box->box distance.
Datum gidx_in(PG_FUNCTION_ARGS)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
static bool gserialized_gist_picksplit_badratios(int *pos, int dims)
Datum gserialized_gist_same(PG_FUNCTION_ARGS)
static int gserialized_gist_picksplit_badratio(int x, int y)
static void gserialized_gist_picksplit_fallback(GistEntryVector *entryvec, GIST_SPLITVEC *v)
static double gidx_distance_leaf_centroid(const GIDX *a, const GIDX *b)
Calculate the centroid->centroid distance between the boxes.
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition: lwpoint.c:57
static float gidx_inter_volume(GIDX *a, GIDX *b)
tuple x
Definition: pixval.py:53
Datum distance(PG_FUNCTION_ARGS)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:102
bool(* gidx_predicate)(GIDX *a, GIDX *b)
static void gserialized_gist_picksplit_addlist(OffsetNumber *list, GIDX **box_union, GIDX *box_current, int *pos, int num)
Datum gserialized_gist_picksplit(PG_FUNCTION_ARGS)
static bool gidx_overlaps(GIDX *a, GIDX *b)
#define FALSE
Definition: dbfopen.c:168
Datum gserialized_overlaps(PG_FUNCTION_ARGS)
double lwgeom_length(const LWGEOM *geom)
Definition: lwgeom.c:1659
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
#define FPeq(A, B)
Definition: box2d.c:37
LWPOINT * lwline_get_lwpoint(const LWLINE *line, int where)
Returns freshly allocated LWPOINT that corresponds to the index where.
Definition: lwline.c:324
LWGEOM * lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:42
Datum gserialized_gist_consistent(PG_FUNCTION_ARGS)
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:856
#define ABOVE(d)
#define TRUE
Definition: dbfopen.c:169
static bool gserialized_gist_consistent_internal(GIDX *key, GIDX *query, StrategyNumber strategy)
Datum gidx_out(PG_FUNCTION_ARGS)
tuple y
Definition: pixval.py:54
This library is the generic geometry handling section of PostGIS.