pgLatLon
view latlon-v0006.c @ 33:ae622355c4d4
Write Earth with capital E in README
| author | jbe |
|---|---|
| date | Sat Sep 24 20:28:28 2016 +0200 (2016-09-24) |
| parents | f765aed73421 |
| children |
line source
2 /*-------------*
3 * C prelude *
4 *-------------*/
6 #include "postgres.h"
7 #include "fmgr.h"
8 #include "libpq/pqformat.h"
9 #include "access/gist.h"
10 #include "access/stratnum.h"
11 #include "utils/array.h"
12 #include <math.h>
14 #ifdef PG_MODULE_MAGIC
15 PG_MODULE_MAGIC;
16 #endif
18 #if INT_MAX < 2147483647
19 #error Expected int type to be at least 32 bit wide
20 #endif
23 /*---------------------------------*
24 * distance calculation on earth *
25 * (using WGS-84 spheroid) *
26 *---------------------------------*/
28 /* WGS-84 spheroid with following parameters:
29 semi-major axis a = 6378137
30 semi-minor axis b = a * (1 - 1/298.257223563)
31 estimated diameter = 2 * (2*a+b)/3
32 */
33 #define PGL_SPHEROID_A 6378137.0 /* semi major axis */
34 #define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */
35 #define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
36 #define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
37 PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
38 ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
39 #define PGL_SUBEPS2 (1.0-PGL_EPS2)
40 #define PGL_DIAMETER ((4.0*PGL_SPHEROID_A + 2.0*PGL_SPHEROID_B) / 3.0)
41 #define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */
42 #define PGL_FADELIMIT (PGL_DIAMETER * M_PI / 6.0) /* 1/6 circumference */
43 #define PGL_MAXDIST (PGL_DIAMETER * M_PI / 2.0) /* maximum distance */
45 /* calculate distance between two points on earth (given in degrees) */
46 static inline double pgl_distance(
47 double lat1, double lon1, double lat2, double lon2
48 ) {
49 float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
50 float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
51 /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
52 /* lon1 = 0 (not used anymore) */
53 lon2 = fabs(lon2-lon1);
54 /* convert to radians (first divide, then multiply) */
55 lat1 = (lat1 / 180.0) * M_PI;
56 lat2 = (lat2 / 180.0) * M_PI;
57 lon2 = (lon2 / 180.0) * M_PI;
58 /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
59 operations (lon2 >= lon1 is already ensured in a previous step) */
60 if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
61 /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
62 of 1.0 */
63 lat1cos = cos(lat1); lat1sin = sin(lat1);
64 lat2cos = cos(lat2); lat2sin = sin(lat2);
65 lon2cos = cos(lon2); lon2sin = sin(lon2);
66 nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
67 nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
68 x1 = nphi1 * lat1cos;
69 z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
70 x2 = nphi2 * lat2cos * lon2cos;
71 y2 = nphi2 * lat2cos * lon2sin;
72 z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
73 /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
74 g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
75 /* convert tunnel distance through scaled ellipsoid to approximated surface
76 distance on original ellipsoid */
77 if (g > 1.0) g = 1.0;
78 s = PGL_DIAMETER * asin(g);
79 /* return result only if small enough to be precise (less than 1/3 of
80 maximum possible distance) */
81 if (s <= PGL_FADELIMIT) return s;
82 /* calculate tunnel distance to antipodal point through scaled ellipsoid */
83 g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
84 /* convert tunnel distance to antipodal point through scaled ellipsoid to
85 approximated surface distance to antipodal point on original ellipsoid */
86 if (g > 1.0) g = 1.0;
87 t = PGL_DIAMETER * asin(g);
88 /* surface distance between original points can now be approximated by
89 substracting antipodal distance from maximum possible distance;
90 return result only if small enough (less than 1/3 of maximum possible
91 distance) */
92 if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
93 /* otherwise crossfade direct and antipodal result to ensure monotonicity */
94 return (
95 (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
96 (s + t - 2*PGL_FADELIMIT)
97 );
98 }
100 /* finite distance that can not be reached on earth */
101 #define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
104 /*--------------------------------*
105 * simple geographic data types *
106 *--------------------------------*/
108 /* point on earth given by latitude and longitude in degrees */
109 /* (type "epoint" in SQL) */
110 typedef struct {
111 double lat; /* between -90 and 90 (both inclusive) */
112 double lon; /* between -180 and 180 (both inclusive) */
113 } pgl_point;
115 /* box delimited by two parallels and two meridians (all in degrees) */
116 /* (type "ebox" in SQL) */
117 typedef struct {
118 double lat_min; /* between -90 and 90 (both inclusive) */
119 double lat_max; /* between -90 and 90 (both inclusive) */
120 double lon_min; /* between -180 and 180 (both inclusive) */
121 double lon_max; /* between -180 and 180 (both inclusive) */
122 /* if lat_min > lat_max, then box is empty */
123 /* if lon_min > lon_max, then 180th meridian is crossed */
124 } pgl_box;
126 /* circle on earth surface (for radial searches with fixed radius) */
127 /* (type "ecircle" in SQL) */
128 typedef struct {
129 pgl_point center;
130 double radius; /* positive (including +0 but excluding -0), or -INFINITY */
131 /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
132 zero radius (0) denotes a single point,
133 a finite radius (0 < radius < INFINITY) denotes a filled circle, and
134 a radius of INFINITY is valid and means complete coverage of earth. */
135 } pgl_circle;
138 /*----------------------------------*
139 * geographic "cluster" data type *
140 *----------------------------------*/
142 /* A cluster is a collection of points, paths, outlines, and polygons. If two
143 polygons in a cluster overlap, the area covered by both polygons does not
144 belong to the cluster. This way, a cluster can be used to describe complex
145 shapes like polygons with holes. Outlines are non-filled polygons. Paths are
146 open by default (i.e. the last point in the list is not connected with the
147 first point in the list). Note that each outline or polygon in a cluster
148 must cover a longitude range of less than 180 degrees to avoid ambiguities.
149 Areas which are larger may be split into multiple polygons. */
151 /* maximum number of points in a cluster */
152 /* (limited to avoid integer overflows, e.g. when allocating memory) */
153 #define PGL_CLUSTER_MAXPOINTS 16777216
155 /* types of cluster entries */
156 #define PGL_ENTRY_POINT 1 /* a point */
157 #define PGL_ENTRY_PATH 2 /* a path from first point to last point */
158 #define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */
159 #define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */
161 /* Entries of a cluster are described by two different structs: pgl_newentry
162 and pgl_entry. The first is used only during construction of a cluster, the
163 second is used in all other cases (e.g. when reading clusters from the
164 database, performing operations, etc). */
166 /* entry for new geographic cluster during construction of that cluster */
167 typedef struct {
168 int32_t entrytype;
169 int32_t npoints;
170 pgl_point *points; /* pointer to an array of points (pgl_point) */
171 } pgl_newentry;
173 /* entry of geographic cluster */
174 typedef struct {
175 int32_t entrytype; /* type of entry: point, path, outline, polygon */
176 int32_t npoints; /* number of stored points (set to 1 for point entry) */
177 int32_t offset; /* offset of pgl_point array from cluster base address */
178 /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
179 } pgl_entry;
181 /* geographic cluster which is a collection of points, (open) paths, polygons,
182 and outlines (non-filled polygons) */
183 typedef struct {
184 char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */
185 int32_t nentries; /* number of stored points */
186 pgl_circle bounding; /* bounding circle */
187 /* Note: bounding circle ensures alignment of pgl_cluster for points */
188 pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */
189 } pgl_cluster;
191 /* macro to determine memory alignment of points */
192 /* (needed to store pgl_point array after entries in pgl_cluster) */
193 typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
194 #define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
196 /* macro to extract a pointer to the array of points of a cluster entry */
197 #define PGL_ENTRY_POINTS(cluster, idx) \
198 ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
200 /* convert pgl_newentry array to pgl_cluster */
201 static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
202 int i; /* index of current entry */
203 int npoints = 0; /* number of points in whole cluster */
204 int entry_npoints; /* number of points in current entry */
205 int points_offset = PGL_POINT_ALIGNMENT * (
206 ( offsetof(pgl_cluster, entries) +
207 nentries * sizeof(pgl_entry) +
208 PGL_POINT_ALIGNMENT - 1
209 ) / PGL_POINT_ALIGNMENT
210 ); /* offset of pgl_point array from base address (considering alignment) */
211 pgl_cluster *cluster; /* new cluster to be returned */
212 /* determine total number of points */
213 for (i=0; i<nentries; i++) npoints += entries[i].npoints;
214 /* allocate memory for cluster (including entries and points) */
215 cluster = palloc(points_offset + npoints * sizeof(pgl_point));
216 /* re-count total number of points to determine offset for each entry */
217 npoints = 0;
218 /* copy entries and points */
219 for (i=0; i<nentries; i++) {
220 /* determine number of points in entry */
221 entry_npoints = entries[i].npoints;
222 /* copy entry */
223 cluster->entries[i].entrytype = entries[i].entrytype;
224 cluster->entries[i].npoints = entry_npoints;
225 /* calculate offset (in bytes) of pgl_point array */
226 cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
227 /* copy points */
228 memcpy(
229 PGL_ENTRY_POINTS(cluster, i),
230 entries[i].points,
231 entry_npoints * sizeof(pgl_point)
232 );
233 /* update total number of points processed */
234 npoints += entry_npoints;
235 }
236 /* set number of entries in cluster */
237 cluster->nentries = nentries;
238 /* set PostgreSQL header for variable sized data */
239 SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
240 /* return newly created cluster */
241 return cluster;
242 }
245 /*----------------------------------------*
246 * C functions on geographic data types *
247 *----------------------------------------*/
249 /* round latitude or longitude to 12 digits after decimal point */
250 static inline double pgl_round(double val) {
251 return round(val * 1e12) / 1e12;
252 }
254 /* compare two points */
255 /* (equality when same point on earth is described, otherwise an arbitrary
256 linear order) */
257 static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
258 double lon1, lon2; /* modified longitudes for special cases */
259 /* use latitude as first ordering criterion */
260 if (point1->lat < point2->lat) return -1;
261 if (point1->lat > point2->lat) return 1;
262 /* determine modified longitudes (considering special case of poles and
263 180th meridian which can be described as W180 or E180) */
264 if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
265 else if (point1->lon == 180) lon1 = -180;
266 else lon1 = point1->lon;
267 if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
268 else if (point2->lon == 180) lon2 = -180;
269 else lon2 = point2->lon;
270 /* use (modified) longitude as secondary ordering criterion */
271 if (lon1 < lon2) return -1;
272 if (lon1 > lon2) return 1;
273 /* no difference found, points are equal */
274 return 0;
275 }
277 /* compare two boxes */
278 /* (equality when same box on earth is described, otherwise an arbitrary linear
279 order) */
280 static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
281 /* two empty boxes are equal, and an empty box is always considered "less
282 than" a non-empty box */
283 if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
284 if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
285 if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
286 /* use southern border as first ordering criterion */
287 if (box1->lat_min < box2->lat_min) return -1;
288 if (box1->lat_min > box2->lat_min) return 1;
289 /* use northern border as second ordering criterion */
290 if (box1->lat_max < box2->lat_max) return -1;
291 if (box1->lat_max > box2->lat_max) return 1;
292 /* use western border as third ordering criterion */
293 if (box1->lon_min < box2->lon_min) return -1;
294 if (box1->lon_min > box2->lon_min) return 1;
295 /* use eastern border as fourth ordering criterion */
296 if (box1->lon_max < box2->lon_max) return -1;
297 if (box1->lon_max > box2->lon_max) return 1;
298 /* no difference found, boxes are equal */
299 return 0;
300 }
302 /* compare two circles */
303 /* (equality when same circle on earth is described, otherwise an arbitrary
304 linear order) */
305 static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
306 /* two circles with same infinite radius (positive or negative infinity) are
307 considered equal independently of center point */
308 if (
309 !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
310 circle1->radius == circle2->radius
311 ) return 0;
312 /* use radius as first ordering criterion */
313 if (circle1->radius < circle2->radius) return -1;
314 if (circle1->radius > circle2->radius) return 1;
315 /* use center point as secondary ordering criterion */
316 return pgl_point_cmp(&(circle1->center), &(circle2->center));
317 }
319 /* set box to empty box*/
320 static void pgl_box_set_empty(pgl_box *box) {
321 box->lat_min = INFINITY;
322 box->lat_max = -INFINITY;
323 box->lon_min = 0;
324 box->lon_max = 0;
325 }
327 /* check if point is inside a box */
328 static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
329 return (
330 point->lat >= box->lat_min && point->lat <= box->lat_max && (
331 (box->lon_min > box->lon_max) ? (
332 /* box crosses 180th meridian */
333 point->lon >= box->lon_min || point->lon <= box->lon_max
334 ) : (
335 /* box does not cross the 180th meridian */
336 point->lon >= box->lon_min && point->lon <= box->lon_max
337 )
338 )
339 );
340 }
342 /* check if two boxes overlap */
343 static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
344 return (
345 box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */
346 ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
347 ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
348 (
349 /* check if one and only one box crosses the 180th meridian */
350 ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
351 ((box2->lon_min > box2->lon_max) ? 1 : 0)
352 ) ? (
353 /* exactly one box crosses the 180th meridian */
354 box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
355 box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
356 ) : (
357 /* no box or both boxes cross the 180th meridian */
358 (
359 (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
360 (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
361 ) ||
362 /* handle W180 == E180 */
363 ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
364 ( box2->lon_min == -180 && box1->lon_max == 180 )
365 )
366 )
367 );
368 }
370 /* check unambiguousness of east/west orientation of cluster entries and set
371 bounding circle of cluster */
372 static bool pgl_finalize_cluster(pgl_cluster *cluster) {
373 int i, j; /* i: index of entry, j: index of point in entry */
374 int npoints; /* number of points in entry */
375 int total_npoints = 0; /* total number of points in cluster */
376 pgl_point *points; /* points in entry */
377 int lon_dir; /* first point of entry west (-1) or east (+1) */
378 double lon_break = 0; /* antipodal longitude of first point in entry */
379 double lon_min, lon_max; /* covered longitude range of entry */
380 double value; /* temporary variable */
381 /* reset bounding circle center to empty circle at 0/0 coordinates */
382 cluster->bounding.center.lat = 0;
383 cluster->bounding.center.lon = 0;
384 cluster->bounding.radius = -INFINITY;
385 /* if cluster is not empty */
386 if (cluster->nentries != 0) {
387 /* iterate over all cluster entries and ensure they each cover a longitude
388 range less than 180 degrees */
389 for (i=0; i<cluster->nentries; i++) {
390 /* get properties of entry */
391 npoints = cluster->entries[i].npoints;
392 points = PGL_ENTRY_POINTS(cluster, i);
393 /* get longitude of first point of entry */
394 value = points[0].lon;
395 /* initialize lon_min and lon_max with longitude of first point */
396 lon_min = value;
397 lon_max = value;
398 /* determine east/west orientation of first point and calculate antipodal
399 longitude (Note: rounding required here) */
400 if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
401 else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); }
402 else lon_dir = 0;
403 /* iterate over all other points in entry */
404 for (j=1; j<npoints; j++) {
405 /* consider longitude wrap-around */
406 value = points[j].lon;
407 if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
408 else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
409 /* update lon_min and lon_max */
410 if (value < lon_min) lon_min = value;
411 else if (value > lon_max) lon_max = value;
412 /* return false if 180 degrees or more are covered */
413 if (lon_max - lon_min >= 180) return false;
414 }
415 }
416 /* iterate over all points of all entries and calculate arbitrary center
417 point for bounding circle (best if center point minimizes the radius,
418 but some error is allowed here) */
419 for (i=0; i<cluster->nentries; i++) {
420 /* get properties of entry */
421 npoints = cluster->entries[i].npoints;
422 points = PGL_ENTRY_POINTS(cluster, i);
423 /* check if first entry */
424 if (i==0) {
425 /* get longitude of first point of first entry in whole cluster */
426 value = points[0].lon;
427 /* initialize lon_min and lon_max with longitude of first point of
428 first entry in whole cluster (used to determine if whole cluster
429 covers a longitude range of 180 degrees or more) */
430 lon_min = value;
431 lon_max = value;
432 /* determine east/west orientation of first point and calculate
433 antipodal longitude (Note: rounding not necessary here) */
434 if (value < 0) { lon_dir = -1; lon_break = value + 180; }
435 else if (value > 0) { lon_dir = 1; lon_break = value - 180; }
436 else lon_dir = 0;
437 }
438 /* iterate over all points in entry */
439 for (j=0; j<npoints; j++) {
440 /* longitude wrap-around (Note: rounding not necessary here) */
441 value = points[j].lon;
442 if (lon_dir < 0 && value > lon_break) value -= 360;
443 else if (lon_dir > 0 && value < lon_break) value += 360;
444 if (value < lon_min) lon_min = value;
445 else if (value > lon_max) lon_max = value;
446 /* set bounding circle to cover whole earth if more than 180 degrees
447 are covered */
448 if (lon_max - lon_min >= 180) {
449 cluster->bounding.center.lat = 0;
450 cluster->bounding.center.lon = 0;
451 cluster->bounding.radius = INFINITY;
452 return true;
453 }
454 /* add point to bounding circle center (for average calculation) */
455 cluster->bounding.center.lat += points[j].lat;
456 cluster->bounding.center.lon += value;
457 }
458 /* count total number of points */
459 total_npoints += npoints;
460 }
461 /* determine average latitude and longitude of cluster */
462 cluster->bounding.center.lat /= total_npoints;
463 cluster->bounding.center.lon /= total_npoints;
464 /* normalize longitude of center of cluster bounding circle */
465 if (cluster->bounding.center.lon < -180) {
466 cluster->bounding.center.lon += 360;
467 }
468 else if (cluster->bounding.center.lon > 180) {
469 cluster->bounding.center.lon -= 360;
470 }
471 /* round bounding circle center (useful if it is used by other functions) */
472 cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
473 cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
474 /* calculate radius of bounding circle */
475 for (i=0; i<cluster->nentries; i++) {
476 npoints = cluster->entries[i].npoints;
477 points = PGL_ENTRY_POINTS(cluster, i);
478 for (j=0; j<npoints; j++) {
479 value = pgl_distance(
480 cluster->bounding.center.lat, cluster->bounding.center.lon,
481 points[j].lat, points[j].lon
482 );
483 if (value > cluster->bounding.radius) cluster->bounding.radius = value;
484 }
485 }
486 }
487 /* return true (east/west orientation is unambiguous) */
488 return true;
489 }
491 /* check if point is inside cluster */
492 /* (if point is on perimeter, then true is returned if and only if
493 strict == false) */
494 static bool pgl_point_in_cluster(
495 pgl_point *point,
496 pgl_cluster *cluster,
497 bool strict
498 ) {
499 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
500 int entrytype; /* type of entry */
501 int npoints; /* number of points in entry */
502 pgl_point *points; /* array of points in entry */
503 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
504 double lon_break = 0; /* antipodal longitude of first vertex */
505 double lat0 = point->lat; /* latitude of point */
506 double lon0; /* (adjusted) longitude of point */
507 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
508 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
509 double lon; /* longitude of intersection */
510 int counter = 0; /* counter for intersections east of point */
511 /* iterate over all entries */
512 for (i=0; i<cluster->nentries; i++) {
513 /* get type of entry */
514 entrytype = cluster->entries[i].entrytype;
515 /* skip all entries but polygons if perimeters are excluded */
516 if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
517 /* get points of entry */
518 npoints = cluster->entries[i].npoints;
519 points = PGL_ENTRY_POINTS(cluster, i);
520 /* determine east/west orientation of first point of entry and calculate
521 antipodal longitude */
522 lon_break = points[0].lon;
523 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
524 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
525 else lon_dir = 0;
526 /* get longitude of point */
527 lon0 = point->lon;
528 /* consider longitude wrap-around for point */
529 if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
530 else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
531 /* iterate over all edges and vertices */
532 for (j=0; j<npoints; j++) {
533 /* return if point is on vertex of polygon */
534 if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
535 /* calculate index of next vertex */
536 k = (j+1) % npoints;
537 /* skip last edge unless entry is (closed) outline or polygon */
538 if (
539 k == 0 &&
540 entrytype != PGL_ENTRY_OUTLINE &&
541 entrytype != PGL_ENTRY_POLYGON
542 ) continue;
543 /* use previously calculated values for lat1 and lon1 if possible */
544 if (j) {
545 lat1 = lat2;
546 lon1 = lon2;
547 } else {
548 /* otherwise get latitude and longitude values of first vertex */
549 lat1 = points[0].lat;
550 lon1 = points[0].lon;
551 /* and consider longitude wrap-around for first vertex */
552 if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
553 else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
554 }
555 /* get latitude and longitude of next vertex */
556 lat2 = points[k].lat;
557 lon2 = points[k].lon;
558 /* consider longitude wrap-around for next vertex */
559 if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
560 else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
561 /* return if point is on horizontal (west to east) edge of polygon */
562 if (
563 lat0 == lat1 && lat0 == lat2 &&
564 ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
565 ) return !strict;
566 /* check if edge crosses east/west line of point */
567 if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
568 /* calculate longitude of intersection */
569 lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
570 /* return if intersection goes (approximately) through point */
571 if (pgl_round(lon) == lon0) return !strict;
572 /* count intersection if east of point and entry is polygon*/
573 if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
574 }
575 }
576 }
577 /* return true if number of intersections is odd */
578 return counter & 1;
579 }
581 /* check if all points of the second cluster are strictly inside the first
582 cluster */
583 static inline bool pgl_all_cluster_points_strictly_in_cluster(
584 pgl_cluster *outer, pgl_cluster *inner
585 ) {
586 int i, j; /* i: entry, j: point in entry */
587 int npoints; /* number of points in entry */
588 pgl_point *points; /* array of points in entry */
589 /* iterate over all entries of "inner" cluster */
590 for (i=0; i<inner->nentries; i++) {
591 /* get properties of entry */
592 npoints = inner->entries[i].npoints;
593 points = PGL_ENTRY_POINTS(inner, i);
594 /* iterate over all points in entry of "inner" cluster */
595 for (j=0; j<npoints; j++) {
596 /* return false if one point of inner cluster is not in outer cluster */
597 if (!pgl_point_in_cluster(points+j, outer, true)) return false;
598 }
599 }
600 /* otherwise return true */
601 return true;
602 }
604 /* check if any point the second cluster is inside the first cluster */
605 static inline bool pgl_any_cluster_points_in_cluster(
606 pgl_cluster *outer, pgl_cluster *inner
607 ) {
608 int i, j; /* i: entry, j: point in entry */
609 int npoints; /* number of points in entry */
610 pgl_point *points; /* array of points in entry */
611 /* iterate over all entries of "inner" cluster */
612 for (i=0; i<inner->nentries; i++) {
613 /* get properties of entry */
614 npoints = inner->entries[i].npoints;
615 points = PGL_ENTRY_POINTS(inner, i);
616 /* iterate over all points in entry of "inner" cluster */
617 for (j=0; j<npoints; j++) {
618 /* return true if one point of inner cluster is in outer cluster */
619 if (pgl_point_in_cluster(points+j, outer, false)) return true;
620 }
621 }
622 /* otherwise return false */
623 return false;
624 }
626 /* check if line segment strictly crosses line (not just touching) */
627 static inline bool pgl_lseg_crosses_line(
628 double seg_x1, double seg_y1, double seg_x2, double seg_y2,
629 double line_x1, double line_y1, double line_x2, double line_y2
630 ) {
631 return (
632 (
633 (seg_x1-line_x1) * (line_y2-line_y1) -
634 (seg_y1-line_y1) * (line_x2-line_x1)
635 ) * (
636 (seg_x2-line_x1) * (line_y2-line_y1) -
637 (seg_y2-line_y1) * (line_x2-line_x1)
638 )
639 ) < 0;
640 }
642 /* check if paths and outlines of two clusters strictly overlap (not just
643 touching) */
644 static bool pgl_outlines_overlap(
645 pgl_cluster *cluster1, pgl_cluster *cluster2
646 ) {
647 int i1, j1, k1; /* i: entry, j: point in entry, k: next point in entry */
648 int i2, j2, k2;
649 int entrytype1, entrytype2; /* type of entry */
650 int npoints1, npoints2; /* number of points in entry */
651 pgl_point *points1; /* array of points in entry of cluster1 */
652 pgl_point *points2; /* array of points in entry of cluster2 */
653 int lon_dir1, lon_dir2; /* first vertex west (-1) or east (+1) */
654 double lon_break1, lon_break2; /* antipodal longitude of first vertex */
655 double lat11, lon11; /* latitude and (adjusted) longitude of vertex */
656 double lat12, lon12; /* latitude and (adjusted) longitude of next vertex */
657 double lat21, lon21; /* latitude and (adjusted) longitudes for cluster2 */
658 double lat22, lon22;
659 double wrapvalue; /* temporary helper value to adjust wrap-around */
660 /* iterate over all entries of cluster1 */
661 for (i1=0; i1<cluster1->nentries; i1++) {
662 /* get properties of entry in cluster1 and skip points */
663 npoints1 = cluster1->entries[i1].npoints;
664 if (npoints1 < 2) continue;
665 entrytype1 = cluster1->entries[i1].entrytype;
666 points1 = PGL_ENTRY_POINTS(cluster1, i1);
667 /* determine east/west orientation of first point and calculate antipodal
668 longitude */
669 lon_break1 = points1[0].lon;
670 if (lon_break1 < 0) {
671 lon_dir1 = -1;
672 lon_break1 = pgl_round(lon_break1 + 180);
673 } else if (lon_break1 > 0) {
674 lon_dir1 = 1;
675 lon_break1 = pgl_round(lon_break1 - 180);
676 } else lon_dir1 = 0;
677 /* iterate over all edges and vertices in cluster1 */
678 for (j1=0; j1<npoints1; j1++) {
679 /* calculate index of next vertex */
680 k1 = (j1+1) % npoints1;
681 /* skip last edge unless entry is (closed) outline or polygon */
682 if (
683 k1 == 0 &&
684 entrytype1 != PGL_ENTRY_OUTLINE &&
685 entrytype1 != PGL_ENTRY_POLYGON
686 ) continue;
687 /* use previously calculated values for lat1 and lon1 if possible */
688 if (j1) {
689 lat11 = lat12;
690 lon11 = lon12;
691 } else {
692 /* otherwise get latitude and longitude values of first vertex */
693 lat11 = points1[0].lat;
694 lon11 = points1[0].lon;
695 /* and consider longitude wrap-around for first vertex */
696 if (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
697 else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
698 }
699 /* get latitude and longitude of next vertex */
700 lat12 = points1[k1].lat;
701 lon12 = points1[k1].lon;
702 /* consider longitude wrap-around for next vertex */
703 if (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
704 else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
705 /* skip degenerated edges */
706 if (lat11 == lat12 && lon11 == lon12) continue;
707 /* iterate over all entries of cluster2 */
708 for (i2=0; i2<cluster2->nentries; i2++) {
709 /* get points and number of points of entry in cluster2 */
710 npoints2 = cluster2->entries[i2].npoints;
711 if (npoints2 < 2) continue;
712 entrytype2 = cluster2->entries[i2].entrytype;
713 points2 = PGL_ENTRY_POINTS(cluster2, i2);
714 /* determine east/west orientation of first point and calculate antipodal
715 longitude */
716 lon_break2 = points2[0].lon;
717 if (lon_break2 < 0) {
718 lon_dir2 = -1;
719 lon_break2 = pgl_round(lon_break2 + 180);
720 } else if (lon_break2 > 0) {
721 lon_dir2 = 1;
722 lon_break2 = pgl_round(lon_break2 - 180);
723 } else lon_dir2 = 0;
724 /* iterate over all edges and vertices in cluster2 */
725 for (j2=0; j2<npoints2; j2++) {
726 /* calculate index of next vertex */
727 k2 = (j2+1) % npoints2;
728 /* skip last edge unless entry is (closed) outline or polygon */
729 if (
730 k2 == 0 &&
731 entrytype2 != PGL_ENTRY_OUTLINE &&
732 entrytype2 != PGL_ENTRY_POLYGON
733 ) continue;
734 /* use previously calculated values for lat1 and lon1 if possible */
735 if (j2) {
736 lat21 = lat22;
737 lon21 = lon22;
738 } else {
739 /* otherwise get latitude and longitude values of first vertex */
740 lat21 = points2[0].lat;
741 lon21 = points2[0].lon;
742 /* and consider longitude wrap-around for first vertex */
743 if (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
744 else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
745 }
746 /* get latitude and longitude of next vertex */
747 lat22 = points2[k2].lat;
748 lon22 = points2[k2].lon;
749 /* consider longitude wrap-around for next vertex */
750 if (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
751 else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
752 /* skip degenerated edges */
753 if (lat21 == lat22 && lon21 == lon22) continue;
754 /* perform another wrap-around where necessary */
755 /* TODO: improve performance of whole wrap-around mechanism */
756 wrapvalue = (lon21 + lon22) - (lon11 + lon12);
757 if (wrapvalue > 360) {
758 lon21 = pgl_round(lon21 - 360);
759 lon22 = pgl_round(lon22 - 360);
760 } else if (wrapvalue < -360) {
761 lon21 = pgl_round(lon21 + 360);
762 lon22 = pgl_round(lon22 + 360);
763 }
764 /* return true if segments overlap */
765 if (
766 pgl_lseg_crosses_line(
767 lat11, lon11, lat12, lon12,
768 lat21, lon21, lat22, lon22
769 ) && pgl_lseg_crosses_line(
770 lat21, lon21, lat22, lon22,
771 lat11, lon11, lat12, lon12
772 )
773 ) {
774 return true;
775 }
776 }
777 }
778 }
779 }
780 /* otherwise return false */
781 return false;
782 }
784 /* check if second cluster is completely contained in first cluster */
785 static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
786 if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
787 if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
788 if (pgl_outlines_overlap(outer, inner)) return false;
789 return true;
790 }
792 /* check if two clusters overlap */
793 static bool pgl_clusters_overlap(
794 pgl_cluster *cluster1, pgl_cluster *cluster2
795 ) {
796 if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
797 if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
798 if (pgl_outlines_overlap(cluster1, cluster2)) return true;
799 return false;
800 }
803 /* calculate (approximate) distance between point and cluster */
804 static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
805 double comp; /* square of compression of meridians */
806 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
807 int entrytype; /* type of entry */
808 int npoints; /* number of points in entry */
809 pgl_point *points; /* array of points in entry */
810 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
811 double lon_break = 0; /* antipodal longitude of first vertex */
812 double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */
813 double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */
814 double lat0 = point->lat; /* latitude of point */
815 double lon0; /* (adjusted) longitude of point */
816 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
817 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
818 double s; /* scalar for vector calculations */
819 double dist; /* distance calculated in one step */
820 double min_dist = INFINITY; /* minimum distance */
821 /* distance is zero if point is contained in cluster */
822 if (pgl_point_in_cluster(point, cluster, false)) return 0;
823 /* calculate (approximate) square compression of meridians */
824 /* TODO: use more exact formula based on WGS-84 */
825 comp = cos((lat0 / 180.0) * M_PI);
826 comp *= comp;
827 /* iterate over all entries */
828 for (i=0; i<cluster->nentries; i++) {
829 /* get properties of entry */
830 entrytype = cluster->entries[i].entrytype;
831 npoints = cluster->entries[i].npoints;
832 points = PGL_ENTRY_POINTS(cluster, i);
833 /* determine east/west orientation of first point of entry and calculate
834 antipodal longitude */
835 lon_break = points[0].lon;
836 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
837 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
838 else lon_dir = 0;
839 /* determine covered longitude range */
840 for (j=0; j<npoints; j++) {
841 /* get longitude of vertex */
842 lon1 = points[j].lon;
843 /* adjust longitude to fix potential wrap-around */
844 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
845 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
846 /* update minimum and maximum longitude of polygon */
847 if (j == 0 || lon1 < lon_min) lon_min = lon1;
848 if (j == 0 || lon1 > lon_max) lon_max = lon1;
849 }
850 /* adjust longitude wrap-around according to full longitude range */
851 lon_break = (lon_max + lon_min) / 2;
852 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
853 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
854 /* get longitude of point */
855 lon0 = point->lon;
856 /* consider longitude wrap-around for point */
857 if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
858 else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
859 /* iterate over all edges and vertices */
860 for (j=0; j<npoints; j++) {
861 /* use previously calculated values for lat1 and lon1 if possible */
862 if (j) {
863 lat1 = lat2;
864 lon1 = lon2;
865 } else {
866 /* otherwise get latitude and longitude values of first vertex */
867 lat1 = points[0].lat;
868 lon1 = points[0].lon;
869 /* and consider longitude wrap-around for first vertex */
870 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
871 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
872 }
873 /* calculate distance to vertex */
874 dist = pgl_distance(lat0, lon0, lat1, lon1);
875 /* store calculated distance if smallest */
876 if (dist < min_dist) min_dist = dist;
877 /* calculate index of next vertex */
878 k = (j+1) % npoints;
879 /* skip last edge unless entry is (closed) outline or polygon */
880 if (
881 k == 0 &&
882 entrytype != PGL_ENTRY_OUTLINE &&
883 entrytype != PGL_ENTRY_POLYGON
884 ) continue;
885 /* get latitude and longitude of next vertex */
886 lat2 = points[k].lat;
887 lon2 = points[k].lon;
888 /* consider longitude wrap-around for next vertex */
889 if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
890 else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
891 /* go to next vertex and edge if edge is degenerated */
892 if (lat1 == lat2 && lon1 == lon2) continue;
893 /* otherwise test if point can be projected onto edge of polygon */
894 s = (
895 ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
896 ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
897 );
898 /* go to next vertex and edge if point cannot be projected */
899 if (!(s > 0 && s < 1)) continue;
900 /* calculate distance from original point to projected point */
901 dist = pgl_distance(
902 lat0, lon0,
903 lat1 + s * (lat2-lat1),
904 lon1 + s * (lon2-lon1)
905 );
906 /* store calculated distance if smallest */
907 if (dist < min_dist) min_dist = dist;
908 }
909 }
910 /* return minimum distance */
911 return min_dist;
912 }
914 /* calculate (approximate) distance between two clusters */
915 static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
916 int i, j; /* i: entry, j: point in entry */
917 int npoints; /* number of points in entry */
918 pgl_point *points; /* array of points in entry */
919 double dist; /* distance calculated in one step */
920 double min_dist = INFINITY; /* minimum distance */
921 /* consider distance from each point in one cluster to the whole other */
922 for (i=0; i<cluster1->nentries; i++) {
923 npoints = cluster1->entries[i].npoints;
924 points = PGL_ENTRY_POINTS(cluster1, i);
925 for (j=0; j<npoints; j++) {
926 dist = pgl_point_cluster_distance(points+j, cluster2);
927 if (dist == 0) return dist;
928 if (dist < min_dist) min_dist = dist;
929 }
930 }
931 /* consider distance from each point in other cluster to the first cluster */
932 for (i=0; i<cluster2->nentries; i++) {
933 npoints = cluster2->entries[i].npoints;
934 points = PGL_ENTRY_POINTS(cluster2, i);
935 for (j=0; j<npoints; j++) {
936 dist = pgl_point_cluster_distance(points+j, cluster1);
937 if (dist == 0) return dist;
938 if (dist < min_dist) min_dist = dist;
939 }
940 }
941 return min_dist;
942 }
944 /* estimator function for distance between box and point */
945 /* always returns a smaller value than actually correct or zero */
946 static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
947 double dlon; /* longitude range of box (delta longitude) */
948 double distance; /* return value */
949 /* return infinity if box is empty */
950 if (box->lat_min > box->lat_max) return INFINITY;
951 /* return zero if point is inside box */
952 if (pgl_point_in_box(point, box)) return 0;
953 /* calculate delta longitude */
954 dlon = box->lon_max - box->lon_min;
955 if (dlon < 0) dlon += 360; /* 180th meridian crossed */
956 /* if delta longitude is greater than 150 degrees, perform safe fall-back */
957 if (dlon > 150) return 0;
958 /* calculate lower limit for distance (formula below requires dlon <= 150) */
959 /* TODO: provide better estimation function to improve performance */
960 distance = (
961 (1.0-1e-14) * pgl_distance(
962 point->lat,
963 point->lon,
964 (box->lat_min + box->lat_max) / 2,
965 box->lon_min + dlon/2
966 ) - pgl_distance(
967 box->lat_min, box->lon_min,
968 box->lat_max, box->lon_max
969 )
970 );
971 /* truncate negative results to zero */
972 if (distance <= 0) distance = 0;
973 /* return result */
974 return distance;
975 }
978 /*-------------------------------------------------*
979 * geographic index based on space-filling curve *
980 *-------------------------------------------------*/
982 /* number of bytes used for geographic (center) position in keys */
983 #define PGL_KEY_LATLON_BYTELEN 7
985 /* maximum reference value for logarithmic size of geographic objects */
986 #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */
988 /* pointer to index key (either pgl_pointkey or pgl_areakey) */
989 typedef unsigned char *pgl_keyptr;
991 /* index key for points (objects with zero area) on the spheroid */
992 /* bit 0..55: interspersed bits of latitude and longitude,
993 bit 56..57: always zero,
994 bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
995 typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
997 /* index key for geographic objects on spheroid with area greater than zero */
998 /* bit 0..55: interspersed bits of latitude and longitude of center point,
999 bit 56: always set to 1,
1000 bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
1001 bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
1002 PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
1003 = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
1004 (with interspersed bits = 0 and node depth = 0) for keys which
1005 cover both empty and non-empty objects */
1007 typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
1009 /* helper macros for reading/writing index keys */
1010 #define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN
1011 #define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1)
1012 #define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8)
1013 #define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1)
1014 #define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1)
1015 #define PGL_AREAKEY_TYPEMASK 0x80
1016 #define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
1017 #define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
1018 ( PGL_KEY_LATLONBIT(key1, n) ^ \
1019 PGL_KEY_LATLONBIT(key2, n) )
1020 #define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
1021 PGL_AREAKEY_TYPEMASK)
1022 #define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
1023 (PGL_AREAKEY_TYPEMASK-1))
1024 #define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET])
1025 #define PGL_KEY_OBJSIZE_EMPTY 126
1026 #define PGL_KEY_OBJSIZE_UNIVERSAL 127
1027 #define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \
1028 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
1029 PGL_KEY_OBJSIZE_EMPTY )
1030 #define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
1031 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
1032 PGL_KEY_OBJSIZE_UNIVERSAL )
1034 /* set area key to match empty objects only */
1035 static void pgl_key_set_empty(pgl_keyptr key) {
1036 memset(key, 0, sizeof(pgl_areakey));
1037 /* Note: setting node depth to maximum is required for picksplit function */
1038 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
1039 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
1040 }
1042 /* set area key to match any object (including empty objects) */
1043 static void pgl_key_set_universal(pgl_keyptr key) {
1044 memset(key, 0, sizeof(pgl_areakey));
1045 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
1046 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
1047 }
1049 /* convert a point on earth into a max-depth key to be used in index */
1050 static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
1051 double lat = point->lat;
1052 double lon = point->lon;
1053 int i;
1054 /* clear latitude and longitude bits */
1055 memset(key, 0, PGL_KEY_LATLON_BYTELEN);
1056 /* set node depth to maximum and type bit to zero */
1057 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
1058 /* iterate over all latitude/longitude bit pairs */
1059 for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
1060 /* determine latitude bit */
1061 if (lat >= 0) {
1062 key[i/4] |= 0x80 >> (2*(i%4));
1063 lat *= 2; lat -= 90;
1064 } else {
1065 lat *= 2; lat += 90;
1066 }
1067 /* determine longitude bit */
1068 if (lon >= 0) {
1069 key[i/4] |= 0x80 >> (2*(i%4)+1);
1070 lon *= 2; lon -= 180;
1071 } else {
1072 lon *= 2; lon += 180;
1073 }
1074 }
1075 }
1077 /* convert a circle on earth into a max-depth key to be used in an index */
1078 static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
1079 /* handle special case of empty circle */
1080 if (circle->radius < 0) {
1081 pgl_key_set_empty(key);
1082 return;
1083 }
1084 /* perform same action as for point keys */
1085 pgl_point_to_key(&(circle->center), key);
1086 /* but overwrite type and node depth to fit area index key */
1087 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
1088 /* check if radius is greater than (or equal to) reference size */
1089 /* (treat equal values as greater values for numerical safety) */
1090 if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
1091 /* if yes, set logarithmic size to zero */
1092 key[PGL_KEY_OBJSIZE_OFFSET] = 0;
1093 } else {
1094 /* otherwise, determine logarithmic size iteratively */
1095 /* (one step is equivalent to a factor of sqrt(2)) */
1096 double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
1097 int objsize = 1;
1098 while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
1099 /* stop when radius is greater than (or equal to) adjusted reference */
1100 /* (treat equal values as greater values for numerical safety) */
1101 if (circle->radius >= reference) break;
1102 reference /= M_SQRT2;
1103 objsize++;
1104 }
1105 /* set logarithmic size to determined value */
1106 key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
1107 }
1108 }
1110 /* check if one key is subkey of another key or vice versa */
1111 static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
1112 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
1113 /* determine smallest depth */
1114 int depth1 = PGL_KEY_NODEDEPTH(key1);
1115 int depth2 = PGL_KEY_NODEDEPTH(key2);
1116 int depth = (depth1 < depth2) ? depth1 : depth2;
1117 /* check if keys are area keys (assuming that both keys have same type) */
1118 if (PGL_KEY_IS_AREAKEY(key1)) {
1119 int j = 0; /* bit offset for logarithmic object size bits */
1120 int k = 0; /* bit offset for latitude and longitude */
1121 /* fetch logarithmic object size information */
1122 int objsize1 = PGL_KEY_OBJSIZE(key1);
1123 int objsize2 = PGL_KEY_OBJSIZE(key2);
1124 /* handle special cases for empty objects (universal and empty keys) */
1125 if (
1126 objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
1127 objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
1128 ) return true;
1129 if (
1130 objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
1131 objsize2 == PGL_KEY_OBJSIZE_EMPTY
1132 ) return objsize1 == objsize2;
1133 /* iterate through key bits */
1134 for (i=0; i<depth; i++) {
1135 /* every second bit is a bit describing the object size */
1136 if (i%2 == 0) {
1137 /* check if object size bit is different in both keys (objsize1 and
1138 objsize2 describe the minimum index when object size bit is set) */
1139 if (
1140 (objsize1 <= j && objsize2 > j) ||
1141 (objsize2 <= j && objsize1 > j)
1142 ) {
1143 /* bit differs, therefore keys are in separate branches */
1144 return false;
1145 }
1146 /* increase bit counter for object size bits */
1147 j++;
1148 }
1149 /* all other bits describe latitude and longitude */
1150 else {
1151 /* check if bit differs in both keys */
1152 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
1153 /* bit differs, therefore keys are in separate branches */
1154 return false;
1155 }
1156 /* increase bit counter for latitude/longitude bits */
1157 k++;
1158 }
1159 }
1160 }
1161 /* if not, keys are point keys */
1162 else {
1163 /* iterate through key bits */
1164 for (i=0; i<depth; i++) {
1165 /* check if bit differs in both keys */
1166 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
1167 /* bit differs, therefore keys are in separate branches */
1168 return false;
1169 }
1170 }
1171 }
1172 /* return true because keys are in the same branch */
1173 return true;
1174 }
1176 /* combine two keys into new key which covers both original keys */
1177 /* (result stored in first argument) */
1178 static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
1179 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
1180 /* determine smallest depth */
1181 int depth1 = PGL_KEY_NODEDEPTH(dst);
1182 int depth2 = PGL_KEY_NODEDEPTH(src);
1183 int depth = (depth1 < depth2) ? depth1 : depth2;
1184 /* check if keys are area keys (assuming that both keys have same type) */
1185 if (PGL_KEY_IS_AREAKEY(dst)) {
1186 pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */
1187 int j = 0; /* bit offset for logarithmic object size bits */
1188 int k = 0; /* bit offset for latitude and longitude */
1189 /* fetch logarithmic object size information */
1190 int objsize1 = PGL_KEY_OBJSIZE(dst);
1191 int objsize2 = PGL_KEY_OBJSIZE(src);
1192 /* handle special cases for empty objects (universal and empty keys) */
1193 if (
1194 objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
1195 objsize2 > PGL_AREAKEY_MAXOBJSIZE
1196 ) {
1197 if (
1198 objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
1199 objsize2 == PGL_KEY_OBJSIZE_EMPTY
1200 ) pgl_key_set_empty(dst);
1201 else pgl_key_set_universal(dst);
1202 return;
1203 }
1204 /* iterate through key bits */
1205 for (i=0; i<depth; i++) {
1206 /* every second bit is a bit describing the object size */
1207 if (i%2 == 0) {
1208 /* increase bit counter for object size bits first */
1209 /* (handy when setting objsize variable) */
1210 j++;
1211 /* check if object size bit is set in neither key */
1212 if (objsize1 >= j && objsize2 >= j) {
1213 /* set objsize in destination buffer to indicate that size bit is
1214 unset in destination buffer at the current bit position */
1215 dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
1216 }
1217 /* break if object size bit is set in one key only */
1218 else if (objsize1 >= j || objsize2 >= j) break;
1219 }
1220 /* all other bits describe latitude and longitude */
1221 else {
1222 /* break if bit differs in both keys */
1223 if (PGL_KEY_LATLONBIT(dst, k)) {
1224 if (!PGL_KEY_LATLONBIT(src, k)) break;
1225 /* but set bit in destination buffer if bit is set in both keys */
1226 dstbuf[k/8] |= 0x80 >> (k%8);
1227 } else if (PGL_KEY_LATLONBIT(src, k)) break;
1228 /* increase bit counter for latitude/longitude bits */
1229 k++;
1230 }
1231 }
1232 /* set common node depth and type bit (type bit = 1) */
1233 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
1234 /* copy contents of destination buffer to first key */
1235 memcpy(dst, dstbuf, sizeof(pgl_areakey));
1236 }
1237 /* if not, keys are point keys */
1238 else {
1239 pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */
1240 /* iterate through key bits */
1241 for (i=0; i<depth; i++) {
1242 /* break if bit differs in both keys */
1243 if (PGL_KEY_LATLONBIT(dst, i)) {
1244 if (!PGL_KEY_LATLONBIT(src, i)) break;
1245 /* but set bit in destination buffer if bit is set in both keys */
1246 dstbuf[i/8] |= 0x80 >> (i%8);
1247 } else if (PGL_KEY_LATLONBIT(src, i)) break;
1248 }
1249 /* set common node depth (type bit = 0) */
1250 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
1251 /* copy contents of destination buffer to first key */
1252 memcpy(dst, dstbuf, sizeof(pgl_pointkey));
1253 }
1254 }
1256 /* determine center(!) boundaries and radius estimation of index key */
1257 static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
1258 int i;
1259 /* determine node depth */
1260 int depth = PGL_KEY_NODEDEPTH(key);
1261 /* center point of possible result */
1262 double lat = 0;
1263 double lon = 0;
1264 /* maximum distance of real center point from key center */
1265 double dlat = 90;
1266 double dlon = 180;
1267 /* maximum radius of contained objects */
1268 double radius = 0; /* always return zero for point index keys */
1269 /* check if key is area key */
1270 if (PGL_KEY_IS_AREAKEY(key)) {
1271 /* get logarithmic object size */
1272 int objsize = PGL_KEY_OBJSIZE(key);
1273 /* handle special cases for empty objects (universal and empty keys) */
1274 if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
1275 pgl_box_set_empty(box);
1276 return 0;
1277 } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
1278 box->lat_min = -90;
1279 box->lat_max = 90;
1280 box->lon_min = -180;
1281 box->lon_max = 180;
1282 return 0; /* any value >= 0 would do */
1283 }
1284 /* calculate maximum possible radius of objects covered by the given key */
1285 if (objsize == 0) radius = INFINITY;
1286 else {
1287 radius = PGL_AREAKEY_REFOBJSIZE;
1288 while (--objsize) radius /= M_SQRT2;
1289 }
1290 /* iterate over latitude and longitude bits in key */
1291 /* (every second bit is a latitude or longitude bit) */
1292 for (i=0; i<depth/2; i++) {
1293 /* check if latitude bit */
1294 if (i%2 == 0) {
1295 /* cut latitude dimension in half */
1296 dlat /= 2;
1297 /* increase center latitude if bit is 1, otherwise decrease */
1298 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1299 else lat -= dlat;
1300 }
1301 /* otherwise longitude bit */
1302 else {
1303 /* cut longitude dimension in half */
1304 dlon /= 2;
1305 /* increase center longitude if bit is 1, otherwise decrease */
1306 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1307 else lon -= dlon;
1308 }
1309 }
1310 }
1311 /* if not, keys are point keys */
1312 else {
1313 /* iterate over all bits in key */
1314 for (i=0; i<depth; i++) {
1315 /* check if latitude bit */
1316 if (i%2 == 0) {
1317 /* cut latitude dimension in half */
1318 dlat /= 2;
1319 /* increase center latitude if bit is 1, otherwise decrease */
1320 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1321 else lat -= dlat;
1322 }
1323 /* otherwise longitude bit */
1324 else {
1325 /* cut longitude dimension in half */
1326 dlon /= 2;
1327 /* increase center longitude if bit is 1, otherwise decrease */
1328 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1329 else lon -= dlon;
1330 }
1331 }
1332 }
1333 /* calculate boundaries from center point and remaining dlat and dlon */
1334 /* (return values through pointer to box) */
1335 box->lat_min = lat - dlat;
1336 box->lat_max = lat + dlat;
1337 box->lon_min = lon - dlon;
1338 box->lon_max = lon + dlon;
1339 /* return radius (as a function return value) */
1340 return radius;
1341 }
1343 /* estimator function for distance between point and index key */
1344 /* always returns a smaller value than actually correct or zero */
1345 static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
1346 pgl_box box; /* center(!) bounding box of area index key */
1347 /* calculate center(!) bounding box and maximum radius of objects covered
1348 by area index key (radius is zero for point index keys) */
1349 double distance = pgl_key_to_box(key, &box);
1350 /* calculate estimated distance between bounding box of center point of
1351 indexed object and point passed as second argument, then substract maximum
1352 radius of objects covered by index key */
1353 distance = pgl_estimate_point_box_distance(point, &box) - distance;
1354 /* truncate negative results to zero */
1355 if (distance <= 0) distance = 0;
1356 /* return result */
1357 return distance;
1358 }
1361 /*---------------------------------*
1362 * helper functions for text I/O *
1363 *---------------------------------*/
1365 #define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */
1367 /* convert floating point number to string (round-trip safe) */
1368 static void pgl_print_float(char *buf, double flt) {
1369 /* check if number is integral */
1370 if (trunc(flt) == flt) {
1371 /* for integral floats use maximum precision */
1372 snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1373 } else {
1374 /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
1375 snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
1376 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
1377 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1378 }
1379 }
1381 /* convert latitude floating point number (in degrees) to string */
1382 static void pgl_print_lat(char *buf, double lat) {
1383 if (signbit(lat)) {
1384 /* treat negative latitudes (including -0) as south */
1385 snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
1386 } else {
1387 /* treat positive latitudes (including +0) as north */
1388 snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
1389 }
1390 }
1392 /* convert longitude floating point number (in degrees) to string */
1393 static void pgl_print_lon(char *buf, double lon) {
1394 if (signbit(lon)) {
1395 /* treat negative longitudes (including -0) as west */
1396 snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
1397 } else {
1398 /* treat positive longitudes (including +0) as east */
1399 snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
1400 }
1401 }
1403 /* bit masks used as return value of pgl_scan() function */
1404 #define PGL_SCAN_NONE 0 /* no value has been parsed */
1405 #define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */
1406 #define PGL_SCAN_LON (1<<1) /* longitude has been parsed */
1407 #define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */
1409 /* parse a coordinate (can be latitude or longitude) */
1410 static int pgl_scan(char **str, double *lat, double *lon) {
1411 double val;
1412 int len;
1413 if (
1414 sscanf(*str, " N %lf %n", &val, &len) ||
1415 sscanf(*str, " n %lf %n", &val, &len)
1416 ) {
1417 *str += len; *lat = val; return PGL_SCAN_LAT;
1418 }
1419 if (
1420 sscanf(*str, " S %lf %n", &val, &len) ||
1421 sscanf(*str, " s %lf %n", &val, &len)
1422 ) {
1423 *str += len; *lat = -val; return PGL_SCAN_LAT;
1424 }
1425 if (
1426 sscanf(*str, " E %lf %n", &val, &len) ||
1427 sscanf(*str, " e %lf %n", &val, &len)
1428 ) {
1429 *str += len; *lon = val; return PGL_SCAN_LON;
1430 }
1431 if (
1432 sscanf(*str, " W %lf %n", &val, &len) ||
1433 sscanf(*str, " w %lf %n", &val, &len)
1434 ) {
1435 *str += len; *lon = -val; return PGL_SCAN_LON;
1436 }
1437 return PGL_SCAN_NONE;
1438 }
1441 /*-----------------*
1442 * SQL functions *
1443 *-----------------*/
1445 /* Note: These function names use "epoint", "ebox", etc. notation here instead
1446 of "point", "box", etc. in order to distinguish them from any previously
1447 defined functions. */
1449 /* function needed for dummy types and/or not implemented features */
1450 PG_FUNCTION_INFO_V1(pgl_notimpl);
1451 Datum pgl_notimpl(PG_FUNCTION_ARGS) {
1452 ereport(ERROR, (errmsg("not implemented by pgLatLon")));
1453 }
1455 /* set point to latitude and longitude (including checks) */
1456 static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
1457 /* reject infinite or NaN values */
1458 if (!isfinite(lat) || !isfinite(lon)) {
1459 ereport(ERROR, (
1460 errcode(ERRCODE_DATA_EXCEPTION),
1461 errmsg("epoint requires finite coordinates")
1462 ));
1463 }
1464 /* check latitude bounds */
1465 if (lat < -90) {
1466 ereport(WARNING, (errmsg("latitude exceeds south pole")));
1467 lat = -90;
1468 } else if (lat > 90) {
1469 ereport(WARNING, (errmsg("latitude exceeds north pole")));
1470 lat = 90;
1471 }
1472 /* check longitude bounds */
1473 if (lon < -180) {
1474 ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
1475 lon += 360 - trunc(lon / 360) * 360;
1476 } else if (lon > 180) {
1477 ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
1478 lon -= 360 + trunc(lon / 360) * 360;
1479 }
1480 /* store rounded latitude/longitude values for round-trip safety */
1481 point->lat = pgl_round(lat);
1482 point->lon = pgl_round(lon);
1483 }
1485 /* create point ("epoint" in SQL) from latitude and longitude */
1486 PG_FUNCTION_INFO_V1(pgl_create_epoint);
1487 Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
1488 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
1489 pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
1490 PG_RETURN_POINTER(point);
1491 }
1493 /* parse point ("epoint" in SQL) */
1494 /* format: '[NS]<float> [EW]<float>' */
1495 PG_FUNCTION_INFO_V1(pgl_epoint_in);
1496 Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
1497 char *str = PG_GETARG_CSTRING(0); /* input string */
1498 char *strptr = str; /* current position within string */
1499 int done = 0; /* bit mask storing if latitude or longitude was read */
1500 double lat, lon; /* parsed values as double precision floats */
1501 pgl_point *point; /* return value (to be palloc'ed) */
1502 /* parse two floats (each latitude or longitude) separated by white-space */
1503 done |= pgl_scan(&strptr, &lat, &lon);
1504 if (strptr != str && isspace(strptr[-1])) {
1505 done |= pgl_scan(&strptr, &lat, &lon);
1506 }
1507 /* require end of string, and latitude and longitude parsed successfully */
1508 if (strptr[0] || done != PGL_SCAN_LATLON) {
1509 ereport(ERROR, (
1510 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1511 errmsg("invalid input syntax for type epoint: \"%s\"", str)
1512 ));
1513 }
1514 /* allocate memory for result */
1515 point = (pgl_point *)palloc(sizeof(pgl_point));
1516 /* set latitude and longitude (and perform checks) */
1517 pgl_epoint_set_latlon(point, lat, lon);
1518 /* return result */
1519 PG_RETURN_POINTER(point);
1520 }
1522 /* create box ("ebox" in SQL) that is empty */
1523 PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
1524 Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
1525 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1526 pgl_box_set_empty(box);
1527 PG_RETURN_POINTER(box);
1528 }
1530 /* set box to given boundaries (including checks) */
1531 static void pgl_ebox_set_boundaries(
1532 pgl_box *box,
1533 double lat_min, double lat_max, double lon_min, double lon_max
1534 ) {
1535 /* if minimum latitude is greater than maximum latitude, return empty box */
1536 if (lat_min > lat_max) {
1537 pgl_box_set_empty(box);
1538 return;
1539 }
1540 /* otherwise reject infinite or NaN values */
1541 if (
1542 !isfinite(lat_min) || !isfinite(lat_max) ||
1543 !isfinite(lon_min) || !isfinite(lon_max)
1544 ) {
1545 ereport(ERROR, (
1546 errcode(ERRCODE_DATA_EXCEPTION),
1547 errmsg("ebox requires finite coordinates")
1548 ));
1549 }
1550 /* check latitude bounds */
1551 if (lat_max < -90) {
1552 ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
1553 lat_max = -90;
1554 } else if (lat_max > 90) {
1555 ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
1556 lat_max = 90;
1557 }
1558 if (lat_min < -90) {
1559 ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
1560 lat_min = -90;
1561 } else if (lat_min > 90) {
1562 ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
1563 lat_min = 90;
1564 }
1565 /* check if all longitudes are included */
1566 if (lon_max - lon_min >= 360) {
1567 if (lon_max - lon_min > 360) ereport(WARNING, (
1568 errmsg("longitude coverage greater than 360 degrees")
1569 ));
1570 lon_min = -180;
1571 lon_max = 180;
1572 } else {
1573 /* normalize longitude bounds */
1574 if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
1575 else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360;
1576 if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
1577 else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360;
1578 }
1579 /* store rounded latitude/longitude values for round-trip safety */
1580 box->lat_min = pgl_round(lat_min);
1581 box->lat_max = pgl_round(lat_max);
1582 box->lon_min = pgl_round(lon_min);
1583 box->lon_max = pgl_round(lon_max);
1584 /* ensure that rounding does not change orientation */
1585 if (lon_min > lon_max && box->lon_min == box->lon_max) {
1586 box->lon_min = -180;
1587 box->lon_max = 180;
1588 }
1589 }
1591 /* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
1592 PG_FUNCTION_INFO_V1(pgl_create_ebox);
1593 Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
1594 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1595 pgl_ebox_set_boundaries(
1596 box,
1597 PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
1598 PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
1599 );
1600 PG_RETURN_POINTER(box);
1601 }
1603 /* create box ("ebox" in SQL) from two points ("epoint"s) */
1604 /* (can not be used to cover a longitude range of more than 120 degrees) */
1605 PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
1606 Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
1607 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
1608 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
1609 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1610 double lat_min, lat_max, lon_min, lon_max;
1611 double dlon; /* longitude range (delta longitude) */
1612 /* order latitude and longitude boundaries */
1613 if (point2->lat < point1->lat) {
1614 lat_min = point2->lat;
1615 lat_max = point1->lat;
1616 } else {
1617 lat_min = point1->lat;
1618 lat_max = point2->lat;
1619 }
1620 if (point2->lon < point1->lon) {
1621 lon_min = point2->lon;
1622 lon_max = point1->lon;
1623 } else {
1624 lon_min = point1->lon;
1625 lon_max = point2->lon;
1626 }
1627 /* calculate longitude range (round to avoid floating point errors) */
1628 dlon = pgl_round(lon_max - lon_min);
1629 /* determine east-west direction */
1630 if (dlon >= 240) {
1631 /* assume that 180th meridian is crossed and swap min/max longitude */
1632 double swap = lon_min; lon_min = lon_max; lon_max = swap;
1633 } else if (dlon > 120) {
1634 /* unclear orientation since delta longitude > 120 */
1635 ereport(ERROR, (
1636 errcode(ERRCODE_DATA_EXCEPTION),
1637 errmsg("can not determine east/west orientation for ebox")
1638 ));
1639 }
1640 /* use boundaries to setup box (and perform checks) */
1641 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1642 /* return result */
1643 PG_RETURN_POINTER(box);
1644 }
1646 /* parse box ("ebox" in SQL) */
1647 /* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
1648 or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
1649 PG_FUNCTION_INFO_V1(pgl_ebox_in);
1650 Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
1651 char *str = PG_GETARG_CSTRING(0); /* input string */
1652 char *str_lower; /* lower case version of input string */
1653 char *strptr; /* current position within string */
1654 int valid; /* number of valid chars */
1655 int done; /* specifies if latitude or longitude was read */
1656 double val; /* temporary variable */
1657 int lat_count = 0; /* count of latitude values parsed */
1658 int lon_count = 0; /* count of longitufde values parsed */
1659 double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */
1660 pgl_box *box; /* return value (to be palloc'ed) */
1661 /* lowercase input */
1662 str_lower = psprintf("%s", str);
1663 for (strptr=str_lower; *strptr; strptr++) {
1664 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1665 }
1666 /* reset reading position to start of (lowercase) string */
1667 strptr = str_lower;
1668 /* check if empty box */
1669 valid = 0;
1670 sscanf(strptr, " empty %n", &valid);
1671 if (valid && strptr[valid] == 0) {
1672 /* allocate and return empty box */
1673 box = (pgl_box *)palloc(sizeof(pgl_box));
1674 pgl_box_set_empty(box);
1675 PG_RETURN_POINTER(box);
1676 }
1677 /* demand four blocks separated by whitespace */
1678 valid = 0;
1679 sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
1680 /* if four blocks separated by whitespace exist, parse those blocks */
1681 if (strptr[valid] == 0) while (strptr[0]) {
1682 /* parse either latitude or longitude (whichever found in input string) */
1683 done = pgl_scan(&strptr, &val, &val);
1684 /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
1685 if (done == PGL_SCAN_LAT) {
1686 if (!lat_count) lat_min = val; else lat_max = val;
1687 lat_count++;
1688 } else if (done == PGL_SCAN_LON) {
1689 if (!lon_count) lon_min = val; else lon_max = val;
1690 lon_count++;
1691 } else {
1692 break;
1693 }
1694 }
1695 /* require end of string, and two latitude and two longitude values */
1696 if (strptr[0] || lat_count != 2 || lon_count != 2) {
1697 ereport(ERROR, (
1698 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1699 errmsg("invalid input syntax for type ebox: \"%s\"", str)
1700 ));
1701 }
1702 /* free lower case string */
1703 pfree(str_lower);
1704 /* order boundaries (maximum greater than minimum) */
1705 if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
1706 if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
1707 /* allocate memory for result */
1708 box = (pgl_box *)palloc(sizeof(pgl_box));
1709 /* set boundaries (and perform checks) */
1710 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1711 /* return result */
1712 PG_RETURN_POINTER(box);
1713 }
1715 /* set circle to given latitude, longitude, and radius (including checks) */
1716 static void pgl_ecircle_set_latlon_radius(
1717 pgl_circle *circle, double lat, double lon, double radius
1718 ) {
1719 /* set center point (including checks) */
1720 pgl_epoint_set_latlon(&(circle->center), lat, lon);
1721 /* handle non-positive radius */
1722 if (isnan(radius)) {
1723 ereport(ERROR, (
1724 errcode(ERRCODE_DATA_EXCEPTION),
1725 errmsg("invalid radius for ecircle")
1726 ));
1727 }
1728 if (radius == 0) radius = 0; /* avoids -0 */
1729 else if (radius < 0) {
1730 if (isfinite(radius)) {
1731 ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
1732 }
1733 radius = -INFINITY;
1734 }
1735 /* store radius (round-trip safety is ensured by pgl_print_float) */
1736 circle->radius = radius;
1737 }
1739 /* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
1740 PG_FUNCTION_INFO_V1(pgl_create_ecircle);
1741 Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
1742 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1743 pgl_ecircle_set_latlon_radius(
1744 circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
1745 );
1746 PG_RETURN_POINTER(circle);
1747 }
1749 /* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
1750 PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
1751 Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
1752 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1753 double radius = PG_GETARG_FLOAT8(1);
1754 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1755 /* set latitude, longitude, radius (and perform checks) */
1756 pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
1757 /* return result */
1758 PG_RETURN_POINTER(circle);
1759 }
1761 /* parse circle ("ecircle" in SQL) */
1762 /* format: '[NS]<float> [EW]<float> <float>' */
1763 PG_FUNCTION_INFO_V1(pgl_ecircle_in);
1764 Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
1765 char *str = PG_GETARG_CSTRING(0); /* input string */
1766 char *strptr = str; /* current position within string */
1767 double lat, lon, radius; /* parsed values as double precision flaots */
1768 int valid = 0; /* number of valid chars */
1769 int done = 0; /* stores if latitude and/or longitude was read */
1770 pgl_circle *circle; /* return value (to be palloc'ed) */
1771 /* demand three blocks separated by whitespace */
1772 sscanf(strptr, " %*s %*s %*s %n", &valid);
1773 /* if three blocks separated by whitespace exist, parse those blocks */
1774 if (strptr[valid] == 0) {
1775 /* parse latitude and longitude */
1776 done |= pgl_scan(&strptr, &lat, &lon);
1777 done |= pgl_scan(&strptr, &lat, &lon);
1778 /* parse radius (while incrementing strptr by number of bytes parsed) */
1779 valid = 0;
1780 if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
1781 }
1782 /* require end of string and both latitude and longitude being parsed */
1783 if (strptr[0] || done != PGL_SCAN_LATLON) {
1784 ereport(ERROR, (
1785 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1786 errmsg("invalid input syntax for type ecircle: \"%s\"", str)
1787 ));
1788 }
1789 /* allocate memory for result */
1790 circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1791 /* set latitude, longitude, radius (and perform checks) */
1792 pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
1793 /* return result */
1794 PG_RETURN_POINTER(circle);
1795 }
1797 /* parse cluster ("ecluster" in SQL) */
1798 PG_FUNCTION_INFO_V1(pgl_ecluster_in);
1799 Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
1800 int i;
1801 char *str = PG_GETARG_CSTRING(0); /* input string */
1802 char *str_lower; /* lower case version of input string */
1803 char *strptr; /* pointer to current reading position of input */
1804 int npoints_total = 0; /* total number of points in cluster */
1805 int nentries = 0; /* total number of entries */
1806 pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */
1807 int entries_buflen = 4; /* maximum number of elements in entries array */
1808 int valid; /* number of valid chars processed */
1809 double lat, lon; /* latitude and longitude of parsed point */
1810 int entrytype; /* current entry type */
1811 int npoints; /* number of points in current entry */
1812 pgl_point *points; /* array of pgl_point for pgl_newentry */
1813 int points_buflen; /* maximum number of elements in points array */
1814 int done; /* return value of pgl_scan function */
1815 pgl_cluster *cluster; /* created cluster */
1816 /* lowercase input */
1817 str_lower = psprintf("%s", str);
1818 for (strptr=str_lower; *strptr; strptr++) {
1819 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1820 }
1821 /* reset reading position to start of (lowercase) string */
1822 strptr = str_lower;
1823 /* allocate initial buffer for entries */
1824 entries = palloc(entries_buflen * sizeof(pgl_newentry));
1825 /* parse until end of string */
1826 while (strptr[0]) {
1827 /* require previous white-space or closing parenthesis before next token */
1828 if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
1829 goto pgl_ecluster_in_error;
1830 }
1831 /* ignore token "empty" */
1832 valid = 0; sscanf(strptr, " empty %n", &valid);
1833 if (valid) { strptr += valid; continue; }
1834 /* test for "point" token */
1835 valid = 0; sscanf(strptr, " point ( %n", &valid);
1836 if (valid) {
1837 strptr += valid;
1838 entrytype = PGL_ENTRY_POINT;
1839 goto pgl_ecluster_in_type_ok;
1840 }
1841 /* test for "path" token */
1842 valid = 0; sscanf(strptr, " path ( %n", &valid);
1843 if (valid) {
1844 strptr += valid;
1845 entrytype = PGL_ENTRY_PATH;
1846 goto pgl_ecluster_in_type_ok;
1847 }
1848 /* test for "outline" token */
1849 valid = 0; sscanf(strptr, " outline ( %n", &valid);
1850 if (valid) {
1851 strptr += valid;
1852 entrytype = PGL_ENTRY_OUTLINE;
1853 goto pgl_ecluster_in_type_ok;
1854 }
1855 /* test for "polygon" token */
1856 valid = 0; sscanf(strptr, " polygon ( %n", &valid);
1857 if (valid) {
1858 strptr += valid;
1859 entrytype = PGL_ENTRY_POLYGON;
1860 goto pgl_ecluster_in_type_ok;
1861 }
1862 /* error if no valid token found */
1863 goto pgl_ecluster_in_error;
1864 pgl_ecluster_in_type_ok:
1865 /* check if pgl_newentry array needs to grow */
1866 if (nentries == entries_buflen) {
1867 pgl_newentry *newbuf;
1868 entries_buflen *= 2;
1869 newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
1870 memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
1871 pfree(entries);
1872 entries = newbuf;
1873 }
1874 /* reset number of points for current entry */
1875 npoints = 0;
1876 /* allocate array for points */
1877 points_buflen = 4;
1878 points = palloc(points_buflen * sizeof(pgl_point));
1879 /* parse until closing parenthesis */
1880 while (strptr[0] != ')') {
1881 /* error on unexpected end of string */
1882 if (strptr[0] == 0) goto pgl_ecluster_in_error;
1883 /* mark neither latitude nor longitude as read */
1884 done = PGL_SCAN_NONE;
1885 /* require white-space before second, third, etc. point */
1886 if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1887 /* scan latitude (or longitude) */
1888 done |= pgl_scan(&strptr, &lat, &lon);
1889 /* require white-space before second coordinate */
1890 if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1891 /* scan longitude (or latitude) */
1892 done |= pgl_scan(&strptr, &lat, &lon);
1893 /* error unless both latitude and longitude were parsed */
1894 if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
1895 /* throw error if number of points is too high */
1896 if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
1897 ereport(ERROR, (
1898 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1899 errmsg(
1900 "too many points for ecluster entry (maximum %i)",
1901 PGL_CLUSTER_MAXPOINTS
1902 )
1903 ));
1904 }
1905 /* check if pgl_point array needs to grow */
1906 if (npoints == points_buflen) {
1907 pgl_point *newbuf;
1908 points_buflen *= 2;
1909 newbuf = palloc(points_buflen * sizeof(pgl_point));
1910 memcpy(newbuf, points, npoints * sizeof(pgl_point));
1911 pfree(points);
1912 points = newbuf;
1913 }
1914 /* append point to pgl_point array (includes checks) */
1915 pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
1916 /* increase total number of points */
1917 npoints_total++;
1918 }
1919 /* error if entry has no points */
1920 if (!npoints) goto pgl_ecluster_in_error;
1921 /* entries with one point are automatically of type "point" */
1922 if (npoints == 1) entrytype = PGL_ENTRY_POINT;
1923 /* if entries have more than one point */
1924 else {
1925 /* throw error if entry type is "point" */
1926 if (entrytype == PGL_ENTRY_POINT) {
1927 ereport(ERROR, (
1928 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1929 errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
1930 ));
1931 }
1932 /* coerce outlines and polygons with more than 2 points to be a path */
1933 if (npoints == 2) entrytype = PGL_ENTRY_PATH;
1934 }
1935 /* append entry to pgl_newentry array */
1936 entries[nentries].entrytype = entrytype;
1937 entries[nentries].npoints = npoints;
1938 entries[nentries].points = points;
1939 nentries++;
1940 /* consume closing parenthesis */
1941 strptr++;
1942 /* consume white-space */
1943 while (isspace(strptr[0])) strptr++;
1944 }
1945 /* free lower case string */
1946 pfree(str_lower);
1947 /* create cluster from pgl_newentry array */
1948 cluster = pgl_new_cluster(nentries, entries);
1949 /* free pgl_newentry array */
1950 for (i=0; i<nentries; i++) pfree(entries[i].points);
1951 pfree(entries);
1952 /* set bounding circle of cluster and check east/west orientation */
1953 if (!pgl_finalize_cluster(cluster)) {
1954 ereport(ERROR, (
1955 errcode(ERRCODE_DATA_EXCEPTION),
1956 errmsg("can not determine east/west orientation for ecluster"),
1957 errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
1958 ));
1959 }
1960 /* return cluster */
1961 PG_RETURN_POINTER(cluster);
1962 /* code to throw error */
1963 pgl_ecluster_in_error:
1964 ereport(ERROR, (
1965 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1966 errmsg("invalid input syntax for type ecluster: \"%s\"", str)
1967 ));
1968 }
1970 /* convert point ("epoint") to string representation */
1971 PG_FUNCTION_INFO_V1(pgl_epoint_out);
1972 Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
1973 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1974 char latstr[PGL_NUMBUFLEN];
1975 char lonstr[PGL_NUMBUFLEN];
1976 pgl_print_lat(latstr, point->lat);
1977 pgl_print_lon(lonstr, point->lon);
1978 PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
1979 }
1981 /* convert box ("ebox") to string representation */
1982 PG_FUNCTION_INFO_V1(pgl_ebox_out);
1983 Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
1984 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
1985 double lon_min = box->lon_min;
1986 double lon_max = box->lon_max;
1987 char lat_min_str[PGL_NUMBUFLEN];
1988 char lat_max_str[PGL_NUMBUFLEN];
1989 char lon_min_str[PGL_NUMBUFLEN];
1990 char lon_max_str[PGL_NUMBUFLEN];
1991 /* return string "empty" if box is set to be empty */
1992 if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
1993 /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
1994 /* (required since pgl_box_in orders the longitude boundaries) */
1995 if (lon_min > lon_max) {
1996 if (lon_min + lon_max >= 0) lon_min -= 360;
1997 else lon_max += 360;
1998 }
1999 /* format and return result */
2000 pgl_print_lat(lat_min_str, box->lat_min);
2001 pgl_print_lat(lat_max_str, box->lat_max);
2002 pgl_print_lon(lon_min_str, lon_min);
2003 pgl_print_lon(lon_max_str, lon_max);
2004 PG_RETURN_CSTRING(psprintf(
2005 "%s %s %s %s",
2006 lat_min_str, lon_min_str, lat_max_str, lon_max_str
2007 ));
2008 }
2010 /* convert circle ("ecircle") to string representation */
2011 PG_FUNCTION_INFO_V1(pgl_ecircle_out);
2012 Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
2013 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2014 char latstr[PGL_NUMBUFLEN];
2015 char lonstr[PGL_NUMBUFLEN];
2016 char radstr[PGL_NUMBUFLEN];
2017 pgl_print_lat(latstr, circle->center.lat);
2018 pgl_print_lon(lonstr, circle->center.lon);
2019 pgl_print_float(radstr, circle->radius);
2020 PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
2021 }
2023 /* convert cluster ("ecluster") to string representation */
2024 PG_FUNCTION_INFO_V1(pgl_ecluster_out);
2025 Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
2026 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2027 char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */
2028 char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */
2029 char ***strings; /* array of array of strings */
2030 char *string; /* string of current token */
2031 char *res, *resptr; /* result and pointer to current write position */
2032 size_t reslen = 1; /* length of result (init with 1 for terminator) */
2033 int npoints; /* number of points of current entry */
2034 int i, j; /* i: entry, j: point in entry */
2035 /* handle empty clusters */
2036 if (cluster->nentries == 0) {
2037 /* free detoasted cluster (if copy) */
2038 PG_FREE_IF_COPY(cluster, 0);
2039 /* return static result */
2040 PG_RETURN_CSTRING("empty");
2041 }
2042 /* allocate array of array of strings */
2043 strings = palloc(cluster->nentries * sizeof(char **));
2044 /* iterate over all entries in cluster */
2045 for (i=0; i<cluster->nentries; i++) {
2046 /* get number of points in entry */
2047 npoints = cluster->entries[i].npoints;
2048 /* allocate array of strings (one string for each point plus two extra) */
2049 strings[i] = palloc((2 + npoints) * sizeof(char *));
2050 /* determine opening string */
2051 switch (cluster->entries[i].entrytype) {
2052 case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break;
2053 case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break;
2054 case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
2055 case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
2056 default: string = (i==0)?"unknown" :" unknown";
2057 }
2058 /* use opening string as first string in array */
2059 strings[i][0] = string;
2060 /* update result length (for allocating result string later) */
2061 reslen += strlen(string);
2062 /* iterate over all points */
2063 for (j=0; j<npoints; j++) {
2064 /* create string representation of point */
2065 pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
2066 pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
2067 string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
2068 /* copy string pointer to string array */
2069 strings[i][j+1] = string;
2070 /* update result length (for allocating result string later) */
2071 reslen += strlen(string);
2072 }
2073 /* use closing parenthesis as last string in array */
2074 strings[i][npoints+1] = ")";
2075 /* update result length (for allocating result string later) */
2076 reslen++;
2077 }
2078 /* allocate result string */
2079 res = palloc(reslen);
2080 /* set write pointer to begin of result string */
2081 resptr = res;
2082 /* copy strings into result string */
2083 for (i=0; i<cluster->nentries; i++) {
2084 npoints = cluster->entries[i].npoints;
2085 for (j=0; j<npoints+2; j++) {
2086 string = strings[i][j];
2087 strcpy(resptr, string);
2088 resptr += strlen(string);
2089 /* free strings allocated by psprintf */
2090 if (j != 0 && j != npoints+1) pfree(string);
2091 }
2092 /* free array of strings */
2093 pfree(strings[i]);
2094 }
2095 /* free array of array of strings */
2096 pfree(strings);
2097 /* free detoasted cluster (if copy) */
2098 PG_FREE_IF_COPY(cluster, 0);
2099 /* return result */
2100 PG_RETURN_CSTRING(res);
2101 }
2103 /* binary input function for point ("epoint") */
2104 PG_FUNCTION_INFO_V1(pgl_epoint_recv);
2105 Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
2106 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2107 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
2108 point->lat = pq_getmsgfloat8(buf);
2109 point->lon = pq_getmsgfloat8(buf);
2110 PG_RETURN_POINTER(point);
2111 }
2113 /* binary input function for box ("ebox") */
2114 PG_FUNCTION_INFO_V1(pgl_ebox_recv);
2115 Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
2116 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2117 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
2118 box->lat_min = pq_getmsgfloat8(buf);
2119 box->lat_max = pq_getmsgfloat8(buf);
2120 box->lon_min = pq_getmsgfloat8(buf);
2121 box->lon_max = pq_getmsgfloat8(buf);
2122 PG_RETURN_POINTER(box);
2123 }
2125 /* binary input function for circle ("ecircle") */
2126 PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
2127 Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
2128 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2129 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
2130 circle->center.lat = pq_getmsgfloat8(buf);
2131 circle->center.lon = pq_getmsgfloat8(buf);
2132 circle->radius = pq_getmsgfloat8(buf);
2133 PG_RETURN_POINTER(circle);
2134 }
2136 /* TODO: binary receive function for cluster */
2138 /* binary output function for point ("epoint") */
2139 PG_FUNCTION_INFO_V1(pgl_epoint_send);
2140 Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
2141 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2142 StringInfoData buf;
2143 pq_begintypsend(&buf);
2144 pq_sendfloat8(&buf, point->lat);
2145 pq_sendfloat8(&buf, point->lon);
2146 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2147 }
2149 /* binary output function for box ("ebox") */
2150 PG_FUNCTION_INFO_V1(pgl_ebox_send);
2151 Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
2152 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2153 StringInfoData buf;
2154 pq_begintypsend(&buf);
2155 pq_sendfloat8(&buf, box->lat_min);
2156 pq_sendfloat8(&buf, box->lat_max);
2157 pq_sendfloat8(&buf, box->lon_min);
2158 pq_sendfloat8(&buf, box->lon_max);
2159 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2160 }
2162 /* binary output function for circle ("ecircle") */
2163 PG_FUNCTION_INFO_V1(pgl_ecircle_send);
2164 Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
2165 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2166 StringInfoData buf;
2167 pq_begintypsend(&buf);
2168 pq_sendfloat8(&buf, circle->center.lat);
2169 pq_sendfloat8(&buf, circle->center.lon);
2170 pq_sendfloat8(&buf, circle->radius);
2171 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2172 }
2174 /* TODO: binary send functions for cluster */
2176 /* cast point ("epoint") to box ("ebox") */
2177 PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
2178 Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
2179 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2180 pgl_box *box = palloc(sizeof(pgl_box));
2181 box->lat_min = point->lat;
2182 box->lat_max = point->lat;
2183 box->lon_min = point->lon;
2184 box->lon_max = point->lon;
2185 PG_RETURN_POINTER(box);
2186 }
2188 /* cast point ("epoint") to circle ("ecircle") */
2189 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
2190 Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
2191 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2192 pgl_circle *circle = palloc(sizeof(pgl_box));
2193 circle->center = *point;
2194 circle->radius = 0;
2195 PG_RETURN_POINTER(circle);
2196 }
2198 /* cast point ("epoint") to cluster ("ecluster") */
2199 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
2200 Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
2201 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2202 pgl_newentry entry;
2203 entry.entrytype = PGL_ENTRY_POINT;
2204 entry.npoints = 1;
2205 entry.points = point;
2206 PG_RETURN_POINTER(pgl_new_cluster(1, &entry));
2207 }
2209 /* cast box ("ebox") to cluster ("ecluster") */
2210 #define pgl_ebox_to_ecluster_macro(i, a, b) \
2211 entries[i].entrytype = PGL_ENTRY_POLYGON; \
2212 entries[i].npoints = 4; \
2213 entries[i].points = points[i]; \
2214 points[i][0].lat = box->lat_min; \
2215 points[i][0].lon = (a); \
2216 points[i][1].lat = box->lat_min; \
2217 points[i][1].lon = (b); \
2218 points[i][2].lat = box->lat_max; \
2219 points[i][2].lon = (b); \
2220 points[i][3].lat = box->lat_max; \
2221 points[i][3].lon = (a);
2222 PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
2223 Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
2224 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2225 double lon, dlon;
2226 int nentries;
2227 pgl_newentry entries[3];
2228 pgl_point points[3][4];
2229 if (box->lat_min > box->lat_max) {
2230 nentries = 0;
2231 } else if (box->lon_min > box->lon_max) {
2232 if (box->lon_min < 0) {
2233 lon = pgl_round((box->lon_min + 180) / 2.0);
2234 nentries = 3;
2235 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2236 pgl_ebox_to_ecluster_macro(1, lon, 180);
2237 pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
2238 } else if (box->lon_max > 0) {
2239 lon = pgl_round((box->lon_max - 180) / 2.0);
2240 nentries = 3;
2241 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2242 pgl_ebox_to_ecluster_macro(1, -180, lon);
2243 pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
2244 } else {
2245 nentries = 2;
2246 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2247 pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
2248 }
2249 } else {
2250 dlon = pgl_round(box->lon_max - box->lon_min);
2251 if (dlon < 180) {
2252 nentries = 1;
2253 pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
2254 } else {
2255 lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
2256 if (
2257 pgl_round(lon - box->lon_min) < 180 &&
2258 pgl_round(box->lon_max - lon) < 180
2259 ) {
2260 nentries = 2;
2261 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2262 pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
2263 } else {
2264 nentries = 3;
2265 pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
2266 pgl_ebox_to_ecluster_macro(1, -60, 60);
2267 pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
2268 }
2269 }
2270 }
2271 PG_RETURN_POINTER(pgl_new_cluster(nentries, entries));
2272 }
2274 /* extract latitude from point ("epoint") */
2275 PG_FUNCTION_INFO_V1(pgl_epoint_lat);
2276 Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
2277 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
2278 }
2280 /* extract longitude from point ("epoint") */
2281 PG_FUNCTION_INFO_V1(pgl_epoint_lon);
2282 Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
2283 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
2284 }
2286 /* extract minimum latitude from box ("ebox") */
2287 PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
2288 Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
2289 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
2290 }
2292 /* extract maximum latitude from box ("ebox") */
2293 PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
2294 Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
2295 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
2296 }
2298 /* extract minimum longitude from box ("ebox") */
2299 PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
2300 Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
2301 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
2302 }
2304 /* extract maximum longitude from box ("ebox") */
2305 PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
2306 Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
2307 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
2308 }
2310 /* extract center point from circle ("ecircle") */
2311 PG_FUNCTION_INFO_V1(pgl_ecircle_center);
2312 Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
2313 PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
2314 }
2316 /* extract radius from circle ("ecircle") */
2317 PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
2318 Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
2319 PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
2320 }
2322 /* check if point is inside box (overlap operator "&&") in SQL */
2323 PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
2324 Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
2325 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2326 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
2327 PG_RETURN_BOOL(pgl_point_in_box(point, box));
2328 }
2330 /* check if point is inside circle (overlap operator "&&") in SQL */
2331 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
2332 Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
2333 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2334 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2335 PG_RETURN_BOOL(
2336 pgl_distance(
2337 point->lat, point->lon,
2338 circle->center.lat, circle->center.lon
2339 ) <= circle->radius
2340 );
2341 }
2343 /* check if point is inside cluster (overlap operator "&&") in SQL */
2344 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
2345 Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
2346 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2347 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2348 bool retval;
2349 /* points outside bounding circle are always assumed to be non-overlapping
2350 (necessary for consistent table and index scans) */
2351 if (
2352 pgl_distance(
2353 point->lat, point->lon,
2354 cluster->bounding.center.lat, cluster->bounding.center.lon
2355 ) > cluster->bounding.radius
2356 ) retval = false;
2357 else retval = pgl_point_in_cluster(point, cluster, false);
2358 PG_FREE_IF_COPY(cluster, 1);
2359 PG_RETURN_BOOL(retval);
2360 }
2362 /* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
2363 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
2364 Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2365 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2366 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2367 bool retval = pgl_distance(
2368 point->lat, point->lon,
2369 cluster->bounding.center.lat, cluster->bounding.center.lon
2370 ) <= cluster->bounding.radius;
2371 PG_FREE_IF_COPY(cluster, 1);
2372 PG_RETURN_BOOL(retval);
2373 }
2375 /* check if two boxes overlap (overlap operator "&&") in SQL */
2376 PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
2377 Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
2378 pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
2379 pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
2380 PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
2381 }
2383 /* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
2384 PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
2385 Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
2386 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2387 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2388 PG_RETURN_BOOL(
2389 pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
2390 );
2391 }
2393 /* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
2394 PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
2395 Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2396 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2397 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2398 bool retval = pgl_estimate_point_box_distance(
2399 &cluster->bounding.center,
2400 box
2401 ) <= cluster->bounding.radius;
2402 PG_FREE_IF_COPY(cluster, 1);
2403 PG_RETURN_BOOL(retval);
2404 }
2406 /* check if two circles overlap (overlap operator "&&") in SQL */
2407 PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
2408 Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
2409 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2410 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2411 PG_RETURN_BOOL(
2412 pgl_distance(
2413 circle1->center.lat, circle1->center.lon,
2414 circle2->center.lat, circle2->center.lon
2415 ) <= circle1->radius + circle2->radius
2416 );
2417 }
2419 /* check if circle and cluster overlap (overlap operator "&&") in SQL */
2420 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
2421 Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
2422 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2423 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2424 bool retval = (
2425 pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
2426 );
2427 PG_FREE_IF_COPY(cluster, 1);
2428 PG_RETURN_BOOL(retval);
2429 }
2431 /* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
2432 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
2433 Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2434 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2435 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2436 bool retval = pgl_distance(
2437 circle->center.lat, circle->center.lon,
2438 cluster->bounding.center.lat, cluster->bounding.center.lon
2439 ) <= circle->radius + cluster->bounding.radius;
2440 PG_FREE_IF_COPY(cluster, 1);
2441 PG_RETURN_BOOL(retval);
2442 }
2444 /* check if two clusters overlap (overlap operator "&&") in SQL */
2445 PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
2446 Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
2447 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2448 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2449 bool retval;
2450 /* clusters with non-touching bounding circles are always assumed to be
2451 non-overlapping (improves performance and is necessary for consistent
2452 table and index scans) */
2453 if (
2454 pgl_distance(
2455 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
2456 cluster2->bounding.center.lat, cluster2->bounding.center.lon
2457 ) > cluster1->bounding.radius + cluster2->bounding.radius
2458 ) retval = false;
2459 else retval = pgl_clusters_overlap(cluster1, cluster2);
2460 PG_FREE_IF_COPY(cluster1, 0);
2461 PG_FREE_IF_COPY(cluster2, 1);
2462 PG_RETURN_BOOL(retval);
2463 }
2465 /* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
2466 PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
2467 Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2468 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2469 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2470 bool retval = pgl_distance(
2471 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
2472 cluster2->bounding.center.lat, cluster2->bounding.center.lon
2473 ) <= cluster1->bounding.radius + cluster2->bounding.radius;
2474 PG_FREE_IF_COPY(cluster1, 0);
2475 PG_FREE_IF_COPY(cluster2, 1);
2476 PG_RETURN_BOOL(retval);
2477 }
2479 /* check if second cluster is in first cluster (cont. operator "@>) in SQL */
2480 PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
2481 Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
2482 pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2483 pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2484 bool retval;
2485 /* clusters with non-touching bounding circles are always assumed to be
2486 non-overlapping (improves performance and is necessary for consistent
2487 table and index scans) */
2488 if (
2489 pgl_distance(
2490 outer->bounding.center.lat, outer->bounding.center.lon,
2491 inner->bounding.center.lat, inner->bounding.center.lon
2492 ) > outer->bounding.radius + inner->bounding.radius
2493 ) retval = false;
2494 else retval = pgl_cluster_in_cluster(outer, inner);
2495 PG_FREE_IF_COPY(outer, 0);
2496 PG_FREE_IF_COPY(inner, 1);
2497 PG_RETURN_BOOL(retval);
2498 }
2500 /* calculate distance between two points ("<->" operator) in SQL */
2501 PG_FUNCTION_INFO_V1(pgl_epoint_distance);
2502 Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
2503 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
2504 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
2505 PG_RETURN_FLOAT8(pgl_distance(
2506 point1->lat, point1->lon, point2->lat, point2->lon
2507 ));
2508 }
2510 /* calculate point to circle distance ("<->" operator) in SQL */
2511 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
2512 Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
2513 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2514 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2515 double distance = pgl_distance(
2516 point->lat, point->lon, circle->center.lat, circle->center.lon
2517 ) - circle->radius;
2518 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2519 }
2521 /* calculate point to cluster distance ("<->" operator) in SQL */
2522 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
2523 Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
2524 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2525 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2526 double distance = pgl_point_cluster_distance(point, cluster);
2527 PG_FREE_IF_COPY(cluster, 1);
2528 PG_RETURN_FLOAT8(distance);
2529 }
2531 /* calculate distance between two circles ("<->" operator) in SQL */
2532 PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
2533 Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
2534 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2535 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2536 double distance = pgl_distance(
2537 circle1->center.lat, circle1->center.lon,
2538 circle2->center.lat, circle2->center.lon
2539 ) - (circle1->radius + circle2->radius);
2540 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2541 }
2543 /* calculate circle to cluster distance ("<->" operator) in SQL */
2544 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
2545 Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
2546 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2547 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2548 double distance = (
2549 pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
2550 );
2551 PG_FREE_IF_COPY(cluster, 1);
2552 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2553 }
2555 /* calculate distance between two clusters ("<->" operator) in SQL */
2556 PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
2557 Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
2558 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2559 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2560 double retval = pgl_cluster_distance(cluster1, cluster2);
2561 PG_FREE_IF_COPY(cluster1, 0);
2562 PG_FREE_IF_COPY(cluster2, 1);
2563 PG_RETURN_FLOAT8(retval);
2564 }
2567 /*-----------------------------------------------------------*
2568 * B-tree comparison operators and index support functions *
2569 *-----------------------------------------------------------*/
2571 /* macro for a B-tree operator (without detoasting) */
2572 #define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
2573 PG_FUNCTION_INFO_V1(func); \
2574 Datum func(PG_FUNCTION_ARGS) { \
2575 type *a = (type *)PG_GETARG_POINTER(0); \
2576 type *b = (type *)PG_GETARG_POINTER(1); \
2577 PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
2578 }
2580 /* macro for a B-tree comparison function (without detoasting) */
2581 #define PGL_BTREE_CMP(func, type, cmpfunc) \
2582 PG_FUNCTION_INFO_V1(func); \
2583 Datum func(PG_FUNCTION_ARGS) { \
2584 type *a = (type *)PG_GETARG_POINTER(0); \
2585 type *b = (type *)PG_GETARG_POINTER(1); \
2586 PG_RETURN_INT32(cmpfunc(a, b)); \
2587 }
2589 /* macro for a B-tree operator (with detoasting) */
2590 #define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
2591 PG_FUNCTION_INFO_V1(func); \
2592 Datum func(PG_FUNCTION_ARGS) { \
2593 bool res; \
2594 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2595 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2596 res = cmpfunc(a, b) oper 0; \
2597 PG_FREE_IF_COPY(a, 0); \
2598 PG_FREE_IF_COPY(b, 1); \
2599 PG_RETURN_BOOL(res); \
2600 }
2602 /* macro for a B-tree comparison function (with detoasting) */
2603 #define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
2604 PG_FUNCTION_INFO_V1(func); \
2605 Datum func(PG_FUNCTION_ARGS) { \
2606 int32_t res; \
2607 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2608 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2609 res = cmpfunc(a, b); \
2610 PG_FREE_IF_COPY(a, 0); \
2611 PG_FREE_IF_COPY(b, 1); \
2612 PG_RETURN_INT32(res); \
2613 }
2615 /* B-tree operators and comparison function for point */
2616 PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
2617 PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
2618 PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
2619 PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
2620 PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
2621 PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
2622 PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
2624 /* B-tree operators and comparison function for box */
2625 PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
2626 PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
2627 PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
2628 PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
2629 PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
2630 PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
2631 PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
2633 /* B-tree operators and comparison function for circle */
2634 PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
2635 PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
2636 PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
2637 PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
2638 PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
2639 PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
2640 PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
2643 /*--------------------------------*
2644 * GiST index support functions *
2645 *--------------------------------*/
2647 /* GiST "consistent" support function */
2648 PG_FUNCTION_INFO_V1(pgl_gist_consistent);
2649 Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
2650 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2651 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
2652 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
2653 bool *recheck = (bool *)PG_GETARG_POINTER(4);
2654 /* demand recheck because index and query methods are lossy */
2655 *recheck = true;
2656 /* strategy number aliases for different operators using the same strategy */
2657 strategy %= 100;
2658 /* strategy number 11: equality of two points */
2659 if (strategy == 11) {
2660 /* query datum is another point */
2661 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2662 /* convert other point to key */
2663 pgl_pointkey querykey;
2664 pgl_point_to_key(query, querykey);
2665 /* return true if both keys overlap */
2666 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2667 }
2668 /* strategy number 13: equality of two circles */
2669 if (strategy == 13) {
2670 /* query datum is another circle */
2671 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2672 /* convert other circle to key */
2673 pgl_areakey querykey;
2674 pgl_circle_to_key(query, querykey);
2675 /* return true if both keys overlap */
2676 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2677 }
2678 /* for all remaining strategies, keys on empty objects produce no match */
2679 /* (check necessary because query radius may be infinite) */
2680 if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
2681 /* strategy number 21: overlapping with point */
2682 if (strategy == 21) {
2683 /* query datum is a point */
2684 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2685 /* return true if estimated distance (allowed to be smaller than real
2686 distance) between index key and point is zero */
2687 PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
2688 }
2689 /* strategy number 22: (point) overlapping with box */
2690 if (strategy == 22) {
2691 /* query datum is a box */
2692 pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
2693 /* determine bounding box of indexed key */
2694 pgl_box keybox;
2695 pgl_key_to_box(key, &keybox);
2696 /* return true if query box overlaps with bounding box of indexed key */
2697 PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
2698 }
2699 /* strategy number 23: overlapping with circle */
2700 if (strategy == 23) {
2701 /* query datum is a circle */
2702 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2703 /* return true if estimated distance (allowed to be smaller than real
2704 distance) between index key and circle center is smaller than radius */
2705 PG_RETURN_BOOL(
2706 pgl_estimate_key_distance(key, &(query->center)) <= query->radius
2707 );
2708 }
2709 /* strategy number 24: overlapping with cluster */
2710 if (strategy == 24) {
2711 bool retval; /* return value */
2712 /* query datum is a cluster */
2713 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2714 /* return true if estimated distance (allowed to be smaller than real
2715 distance) between index key and circle center is smaller than radius */
2716 retval = (
2717 pgl_estimate_key_distance(key, &(query->bounding.center)) <=
2718 query->bounding.radius
2719 );
2720 PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */
2721 PG_RETURN_BOOL(retval);
2722 }
2723 /* throw error for any unknown strategy number */
2724 elog(ERROR, "unrecognized strategy number: %d", strategy);
2725 }
2727 /* GiST "union" support function */
2728 PG_FUNCTION_INFO_V1(pgl_gist_union);
2729 Datum pgl_gist_union(PG_FUNCTION_ARGS) {
2730 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2731 pgl_keyptr out; /* return value (to be palloc'ed) */
2732 int i;
2733 /* determine key size */
2734 size_t keysize = PGL_KEY_IS_AREAKEY(
2735 (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
2736 ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
2737 /* begin with first key as result */
2738 out = palloc(keysize);
2739 memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
2740 /* unite current result with second, third, etc. key */
2741 for (i=1; i<entryvec->n; i++) {
2742 pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
2743 }
2744 /* return result */
2745 PG_RETURN_POINTER(out);
2746 }
2748 /* GiST "compress" support function for indicis on points */
2749 PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
2750 Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
2751 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2752 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2753 /* only transform new leaves */
2754 if (entry->leafkey) {
2755 /* get point to be transformed */
2756 pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
2757 /* allocate memory for key */
2758 pgl_keyptr key = palloc(sizeof(pgl_pointkey));
2759 /* transform point to key */
2760 pgl_point_to_key(point, key);
2761 /* create new GISTENTRY structure as return value */
2762 retval = palloc(sizeof(GISTENTRY));
2763 gistentryinit(
2764 *retval, PointerGetDatum(key),
2765 entry->rel, entry->page, entry->offset, FALSE
2766 );
2767 } else {
2768 /* inner nodes have already been transformed */
2769 retval = entry;
2770 }
2771 /* return pointer to old or new GISTENTRY structure */
2772 PG_RETURN_POINTER(retval);
2773 }
2775 /* GiST "compress" support function for indicis on circles */
2776 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
2777 Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
2778 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2779 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2780 /* only transform new leaves */
2781 if (entry->leafkey) {
2782 /* get circle to be transformed */
2783 pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
2784 /* allocate memory for key */
2785 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2786 /* transform circle to key */
2787 pgl_circle_to_key(circle, key);
2788 /* create new GISTENTRY structure as return value */
2789 retval = palloc(sizeof(GISTENTRY));
2790 gistentryinit(
2791 *retval, PointerGetDatum(key),
2792 entry->rel, entry->page, entry->offset, FALSE
2793 );
2794 } else {
2795 /* inner nodes have already been transformed */
2796 retval = entry;
2797 }
2798 /* return pointer to old or new GISTENTRY structure */
2799 PG_RETURN_POINTER(retval);
2800 }
2802 /* GiST "compress" support function for indices on clusters */
2803 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
2804 Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
2805 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2806 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2807 /* only transform new leaves */
2808 if (entry->leafkey) {
2809 /* get cluster to be transformed (detoasting necessary!) */
2810 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
2811 /* allocate memory for key */
2812 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2813 /* transform cluster to key */
2814 pgl_circle_to_key(&(cluster->bounding), key);
2815 /* create new GISTENTRY structure as return value */
2816 retval = palloc(sizeof(GISTENTRY));
2817 gistentryinit(
2818 *retval, PointerGetDatum(key),
2819 entry->rel, entry->page, entry->offset, FALSE
2820 );
2821 /* free detoasted datum */
2822 if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
2823 } else {
2824 /* inner nodes have already been transformed */
2825 retval = entry;
2826 }
2827 /* return pointer to old or new GISTENTRY structure */
2828 PG_RETURN_POINTER(retval);
2829 }
2831 /* GiST "decompress" support function for indices */
2832 PG_FUNCTION_INFO_V1(pgl_gist_decompress);
2833 Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
2834 /* return passed pointer without transformation */
2835 PG_RETURN_POINTER(PG_GETARG_POINTER(0));
2836 }
2838 /* GiST "penalty" support function */
2839 PG_FUNCTION_INFO_V1(pgl_gist_penalty);
2840 Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
2841 GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
2842 GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
2843 float *penalty = (float *)PG_GETARG_POINTER(2);
2844 /* get original key and key to insert */
2845 pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
2846 pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
2847 /* copy original key */
2848 union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
2849 if (PGL_KEY_IS_AREAKEY(orig)) {
2850 memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
2851 } else {
2852 memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
2853 }
2854 /* calculate union of both keys */
2855 pgl_unite_keys((pgl_keyptr)&union_key, new);
2856 /* penalty equal to reduction of key length (logarithm of added area) */
2857 /* (return value by setting referenced value and returning pointer) */
2858 *penalty = (
2859 PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
2860 );
2861 PG_RETURN_POINTER(penalty);
2862 }
2864 /* GiST "picksplit" support function */
2865 PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
2866 Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
2867 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2868 GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
2869 OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */
2870 union {
2871 pgl_pointkey pointkey;
2872 pgl_areakey areakey;
2873 } union_all; /* union of all keys (to be calculated from scratch)
2874 (later cut in half) */
2875 int is_areakey = PGL_KEY_IS_AREAKEY(
2876 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
2877 );
2878 int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
2879 pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */
2880 pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */
2881 pgl_keyptr key; /* current key to be processed */
2882 /* allocate memory for array of left and right keys, set counts to zero */
2883 v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2884 v->spl_nleft = 0;
2885 v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2886 v->spl_nright = 0;
2887 /* calculate union of all keys from scratch */
2888 memcpy(
2889 (pgl_keyptr)&union_all,
2890 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
2891 keysize
2892 );
2893 for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
2894 pgl_unite_keys(
2895 (pgl_keyptr)&union_all,
2896 (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
2897 );
2898 }
2899 /* check if trivial split is necessary due to exhausted key length */
2900 /* (Note: keys for empty objects must have node depth set to maximum) */
2901 if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
2902 is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
2903 )) {
2904 /* half of all keys go left */
2905 for (
2906 i=FirstOffsetNumber;
2907 i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2908 i=OffsetNumberNext(i)
2909 ) {
2910 /* pointer to current key */
2911 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2912 /* update unionL */
2913 /* check if key is first key that goes left */
2914 if (!v->spl_nleft) {
2915 /* first key that goes left is just copied to unionL */
2916 memcpy(unionL, key, keysize);
2917 } else {
2918 /* unite current value and next key */
2919 pgl_unite_keys(unionL, key);
2920 }
2921 /* append offset number to list of keys that go left */
2922 v->spl_left[v->spl_nleft++] = i;
2923 }
2924 /* other half goes right */
2925 for (
2926 i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2927 i<entryvec->n;
2928 i=OffsetNumberNext(i)
2929 ) {
2930 /* pointer to current key */
2931 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2932 /* update unionR */
2933 /* check if key is first key that goes right */
2934 if (!v->spl_nright) {
2935 /* first key that goes right is just copied to unionR */
2936 memcpy(unionR, key, keysize);
2937 } else {
2938 /* unite current value and next key */
2939 pgl_unite_keys(unionR, key);
2940 }
2941 /* append offset number to list of keys that go right */
2942 v->spl_right[v->spl_nright++] = i;
2943 }
2944 }
2945 /* otherwise, a non-trivial split is possible */
2946 else {
2947 /* cut covered area in half */
2948 /* (union_all then refers to area of keys that go left) */
2949 /* check if union of all keys covers empty and non-empty objects */
2950 if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
2951 /* if yes, split into empty and non-empty objects */
2952 pgl_key_set_empty((pgl_keyptr)&union_all);
2953 } else {
2954 /* otherwise split by next bit */
2955 ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
2956 /* NOTE: type bit conserved */
2957 }
2958 /* determine for each key if it goes left or right */
2959 for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
2960 /* pointer to current key */
2961 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2962 /* keys within one half of the area go left */
2963 if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
2964 /* update unionL */
2965 /* check if key is first key that goes left */
2966 if (!v->spl_nleft) {
2967 /* first key that goes left is just copied to unionL */
2968 memcpy(unionL, key, keysize);
2969 } else {
2970 /* unite current value of unionL and processed key */
2971 pgl_unite_keys(unionL, key);
2972 }
2973 /* append offset number to list of keys that go left */
2974 v->spl_left[v->spl_nleft++] = i;
2975 }
2976 /* the other keys go right */
2977 else {
2978 /* update unionR */
2979 /* check if key is first key that goes right */
2980 if (!v->spl_nright) {
2981 /* first key that goes right is just copied to unionR */
2982 memcpy(unionR, key, keysize);
2983 } else {
2984 /* unite current value of unionR and processed key */
2985 pgl_unite_keys(unionR, key);
2986 }
2987 /* append offset number to list of keys that go right */
2988 v->spl_right[v->spl_nright++] = i;
2989 }
2990 }
2991 }
2992 /* store unions in return value */
2993 v->spl_ldatum = PointerGetDatum(unionL);
2994 v->spl_rdatum = PointerGetDatum(unionR);
2995 /* return all results */
2996 PG_RETURN_POINTER(v);
2997 }
2999 /* GiST "same"/"equal" support function */
3000 PG_FUNCTION_INFO_V1(pgl_gist_same);
3001 Datum pgl_gist_same(PG_FUNCTION_ARGS) {
3002 pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
3003 pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
3004 bool *boolptr = (bool *)PG_GETARG_POINTER(2);
3005 /* two keys are equal if they are binary equal */
3006 /* (return result by setting referenced boolean and returning pointer) */
3007 *boolptr = !memcmp(
3008 key1,
3009 key2,
3010 PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
3011 );
3012 PG_RETURN_POINTER(boolptr);
3013 }
3015 /* GiST "distance" support function */
3016 PG_FUNCTION_INFO_V1(pgl_gist_distance);
3017 Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
3018 GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
3019 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
3020 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
3021 bool *recheck = (bool *)PG_GETARG_POINTER(4);
3022 double distance; /* return value */
3023 /* demand recheck because distance is just an estimation */
3024 /* (real distance may be bigger) */
3025 *recheck = true;
3026 /* strategy number aliases for different operators using the same strategy */
3027 strategy %= 100;
3028 /* strategy number 31: distance to point */
3029 if (strategy == 31) {
3030 /* query datum is a point */
3031 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
3032 /* use pgl_estimate_pointkey_distance() function to compute result */
3033 distance = pgl_estimate_key_distance(key, query);
3034 /* avoid infinity (reserved!) */
3035 if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3036 /* return result */
3037 PG_RETURN_FLOAT8(distance);
3038 }
3039 /* strategy number 33: distance to circle */
3040 if (strategy == 33) {
3041 /* query datum is a circle */
3042 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
3043 /* estimate distance to circle center and substract circle radius */
3044 distance = (
3045 pgl_estimate_key_distance(key, &(query->center)) - query->radius
3046 );
3047 /* convert non-positive values to zero and avoid infinity (reserved!) */
3048 if (distance <= 0) distance = 0;
3049 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3050 /* return result */
3051 PG_RETURN_FLOAT8(distance);
3052 }
3053 /* strategy number 34: distance to cluster */
3054 if (strategy == 34) {
3055 /* query datum is a cluster */
3056 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3057 /* estimate distance to bounding center and substract bounding radius */
3058 distance = (
3059 pgl_estimate_key_distance(key, &(query->bounding.center)) -
3060 query->bounding.radius
3061 );
3062 /* convert non-positive values to zero and avoid infinity (reserved!) */
3063 if (distance <= 0) distance = 0;
3064 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3065 /* free detoasted cluster (if copy) */
3066 PG_FREE_IF_COPY(query, 1);
3067 /* return result */
3068 PG_RETURN_FLOAT8(distance);
3069 }
3070 /* throw error for any unknown strategy number */
3071 elog(ERROR, "unrecognized strategy number: %d", strategy);
3072 }