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