pgLatLon
view latlon-v0007.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 | bdbec73dc8ff |
| 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 comp = cos((lat0 / 180.0) * M_PI);
825 comp *= comp;
826 /* calculate exact square compression of meridians */
827 comp *= (
828 (1.0 - PGL_EPS2 * (1.0-comp)) *
829 (1.0 - PGL_EPS2 * (1.0-comp)) /
830 (PGL_SUBEPS2 * PGL_SUBEPS2)
831 );
832 /* iterate over all entries */
833 for (i=0; i<cluster->nentries; i++) {
834 /* get properties of entry */
835 entrytype = cluster->entries[i].entrytype;
836 npoints = cluster->entries[i].npoints;
837 points = PGL_ENTRY_POINTS(cluster, i);
838 /* determine east/west orientation of first point of entry and calculate
839 antipodal longitude */
840 lon_break = points[0].lon;
841 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
842 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
843 else lon_dir = 0;
844 /* determine covered longitude range */
845 for (j=0; j<npoints; j++) {
846 /* get longitude of vertex */
847 lon1 = points[j].lon;
848 /* adjust longitude to fix potential wrap-around */
849 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
850 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
851 /* update minimum and maximum longitude of polygon */
852 if (j == 0 || lon1 < lon_min) lon_min = lon1;
853 if (j == 0 || lon1 > lon_max) lon_max = lon1;
854 }
855 /* adjust longitude wrap-around according to full longitude range */
856 lon_break = (lon_max + lon_min) / 2;
857 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
858 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
859 /* get longitude of point */
860 lon0 = point->lon;
861 /* consider longitude wrap-around for point */
862 if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
863 else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
864 /* iterate over all edges and vertices */
865 for (j=0; j<npoints; j++) {
866 /* use previously calculated values for lat1 and lon1 if possible */
867 if (j) {
868 lat1 = lat2;
869 lon1 = lon2;
870 } else {
871 /* otherwise get latitude and longitude values of first vertex */
872 lat1 = points[0].lat;
873 lon1 = points[0].lon;
874 /* and consider longitude wrap-around for first vertex */
875 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
876 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
877 }
878 /* calculate distance to vertex */
879 dist = pgl_distance(lat0, lon0, lat1, lon1);
880 /* store calculated distance if smallest */
881 if (dist < min_dist) min_dist = dist;
882 /* calculate index of next vertex */
883 k = (j+1) % npoints;
884 /* skip last edge unless entry is (closed) outline or polygon */
885 if (
886 k == 0 &&
887 entrytype != PGL_ENTRY_OUTLINE &&
888 entrytype != PGL_ENTRY_POLYGON
889 ) continue;
890 /* get latitude and longitude of next vertex */
891 lat2 = points[k].lat;
892 lon2 = points[k].lon;
893 /* consider longitude wrap-around for next vertex */
894 if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
895 else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
896 /* go to next vertex and edge if edge is degenerated */
897 if (lat1 == lat2 && lon1 == lon2) continue;
898 /* otherwise test if point can be projected onto edge of polygon */
899 s = (
900 ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
901 ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
902 );
903 /* go to next vertex and edge if point cannot be projected */
904 if (!(s > 0 && s < 1)) continue;
905 /* calculate distance from original point to projected point */
906 dist = pgl_distance(
907 lat0, lon0,
908 lat1 + s * (lat2-lat1),
909 lon1 + s * (lon2-lon1)
910 );
911 /* store calculated distance if smallest */
912 if (dist < min_dist) min_dist = dist;
913 }
914 }
915 /* return minimum distance */
916 return min_dist;
917 }
919 /* calculate (approximate) distance between two clusters */
920 static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
921 int i, j; /* i: entry, j: point in entry */
922 int npoints; /* number of points in entry */
923 pgl_point *points; /* array of points in entry */
924 double dist; /* distance calculated in one step */
925 double min_dist = INFINITY; /* minimum distance */
926 /* consider distance from each point in one cluster to the whole other */
927 for (i=0; i<cluster1->nentries; i++) {
928 npoints = cluster1->entries[i].npoints;
929 points = PGL_ENTRY_POINTS(cluster1, i);
930 for (j=0; j<npoints; j++) {
931 dist = pgl_point_cluster_distance(points+j, cluster2);
932 if (dist == 0) return dist;
933 if (dist < min_dist) min_dist = dist;
934 }
935 }
936 /* consider distance from each point in other cluster to the first cluster */
937 for (i=0; i<cluster2->nentries; i++) {
938 npoints = cluster2->entries[i].npoints;
939 points = PGL_ENTRY_POINTS(cluster2, i);
940 for (j=0; j<npoints; j++) {
941 dist = pgl_point_cluster_distance(points+j, cluster1);
942 if (dist == 0) return dist;
943 if (dist < min_dist) min_dist = dist;
944 }
945 }
946 return min_dist;
947 }
949 /* estimator function for distance between box and point */
950 /* always returns a smaller value than actually correct or zero */
951 static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
952 double dlon; /* longitude range of box (delta longitude) */
953 double distance; /* return value */
954 /* return infinity if box is empty */
955 if (box->lat_min > box->lat_max) return INFINITY;
956 /* return zero if point is inside box */
957 if (pgl_point_in_box(point, box)) return 0;
958 /* calculate delta longitude */
959 dlon = box->lon_max - box->lon_min;
960 if (dlon < 0) dlon += 360; /* 180th meridian crossed */
961 /* if delta longitude is greater than 150 degrees, perform safe fall-back */
962 if (dlon > 150) return 0;
963 /* calculate lower limit for distance (formula below requires dlon <= 150) */
964 /* TODO: provide better estimation function to improve performance */
965 distance = (
966 (1.0-1e-14) * pgl_distance(
967 point->lat,
968 point->lon,
969 (box->lat_min + box->lat_max) / 2,
970 box->lon_min + dlon/2
971 ) - pgl_distance(
972 box->lat_min, box->lon_min,
973 box->lat_max, box->lon_max
974 )
975 );
976 /* truncate negative results to zero */
977 if (distance <= 0) distance = 0;
978 /* return result */
979 return distance;
980 }
983 /*-------------------------------------------------*
984 * geographic index based on space-filling curve *
985 *-------------------------------------------------*/
987 /* number of bytes used for geographic (center) position in keys */
988 #define PGL_KEY_LATLON_BYTELEN 7
990 /* maximum reference value for logarithmic size of geographic objects */
991 #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */
993 /* pointer to index key (either pgl_pointkey or pgl_areakey) */
994 typedef unsigned char *pgl_keyptr;
996 /* index key for points (objects with zero area) on the spheroid */
997 /* bit 0..55: interspersed bits of latitude and longitude,
998 bit 56..57: always zero,
999 bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
1000 typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
1002 /* index key for geographic objects on spheroid with area greater than zero */
1003 /* bit 0..55: interspersed bits of latitude and longitude of center point,
1004 bit 56: always set to 1,
1005 bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
1006 bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
1007 PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
1008 = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
1009 (with interspersed bits = 0 and node depth = 0) for keys which
1010 cover both empty and non-empty objects */
1012 typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
1014 /* helper macros for reading/writing index keys */
1015 #define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN
1016 #define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1)
1017 #define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8)
1018 #define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1)
1019 #define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1)
1020 #define PGL_AREAKEY_TYPEMASK 0x80
1021 #define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
1022 #define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
1023 ( PGL_KEY_LATLONBIT(key1, n) ^ \
1024 PGL_KEY_LATLONBIT(key2, n) )
1025 #define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
1026 PGL_AREAKEY_TYPEMASK)
1027 #define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
1028 (PGL_AREAKEY_TYPEMASK-1))
1029 #define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET])
1030 #define PGL_KEY_OBJSIZE_EMPTY 126
1031 #define PGL_KEY_OBJSIZE_UNIVERSAL 127
1032 #define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \
1033 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
1034 PGL_KEY_OBJSIZE_EMPTY )
1035 #define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
1036 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
1037 PGL_KEY_OBJSIZE_UNIVERSAL )
1039 /* set area key to match empty objects only */
1040 static void pgl_key_set_empty(pgl_keyptr key) {
1041 memset(key, 0, sizeof(pgl_areakey));
1042 /* Note: setting node depth to maximum is required for picksplit function */
1043 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
1044 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
1045 }
1047 /* set area key to match any object (including empty objects) */
1048 static void pgl_key_set_universal(pgl_keyptr key) {
1049 memset(key, 0, sizeof(pgl_areakey));
1050 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
1051 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
1052 }
1054 /* convert a point on earth into a max-depth key to be used in index */
1055 static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
1056 double lat = point->lat;
1057 double lon = point->lon;
1058 int i;
1059 /* clear latitude and longitude bits */
1060 memset(key, 0, PGL_KEY_LATLON_BYTELEN);
1061 /* set node depth to maximum and type bit to zero */
1062 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
1063 /* iterate over all latitude/longitude bit pairs */
1064 for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
1065 /* determine latitude bit */
1066 if (lat >= 0) {
1067 key[i/4] |= 0x80 >> (2*(i%4));
1068 lat *= 2; lat -= 90;
1069 } else {
1070 lat *= 2; lat += 90;
1071 }
1072 /* determine longitude bit */
1073 if (lon >= 0) {
1074 key[i/4] |= 0x80 >> (2*(i%4)+1);
1075 lon *= 2; lon -= 180;
1076 } else {
1077 lon *= 2; lon += 180;
1078 }
1079 }
1080 }
1082 /* convert a circle on earth into a max-depth key to be used in an index */
1083 static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
1084 /* handle special case of empty circle */
1085 if (circle->radius < 0) {
1086 pgl_key_set_empty(key);
1087 return;
1088 }
1089 /* perform same action as for point keys */
1090 pgl_point_to_key(&(circle->center), key);
1091 /* but overwrite type and node depth to fit area index key */
1092 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
1093 /* check if radius is greater than (or equal to) reference size */
1094 /* (treat equal values as greater values for numerical safety) */
1095 if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
1096 /* if yes, set logarithmic size to zero */
1097 key[PGL_KEY_OBJSIZE_OFFSET] = 0;
1098 } else {
1099 /* otherwise, determine logarithmic size iteratively */
1100 /* (one step is equivalent to a factor of sqrt(2)) */
1101 double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
1102 int objsize = 1;
1103 while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
1104 /* stop when radius is greater than (or equal to) adjusted reference */
1105 /* (treat equal values as greater values for numerical safety) */
1106 if (circle->radius >= reference) break;
1107 reference /= M_SQRT2;
1108 objsize++;
1109 }
1110 /* set logarithmic size to determined value */
1111 key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
1112 }
1113 }
1115 /* check if one key is subkey of another key or vice versa */
1116 static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
1117 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
1118 /* determine smallest depth */
1119 int depth1 = PGL_KEY_NODEDEPTH(key1);
1120 int depth2 = PGL_KEY_NODEDEPTH(key2);
1121 int depth = (depth1 < depth2) ? depth1 : depth2;
1122 /* check if keys are area keys (assuming that both keys have same type) */
1123 if (PGL_KEY_IS_AREAKEY(key1)) {
1124 int j = 0; /* bit offset for logarithmic object size bits */
1125 int k = 0; /* bit offset for latitude and longitude */
1126 /* fetch logarithmic object size information */
1127 int objsize1 = PGL_KEY_OBJSIZE(key1);
1128 int objsize2 = PGL_KEY_OBJSIZE(key2);
1129 /* handle special cases for empty objects (universal and empty keys) */
1130 if (
1131 objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
1132 objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
1133 ) return true;
1134 if (
1135 objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
1136 objsize2 == PGL_KEY_OBJSIZE_EMPTY
1137 ) return objsize1 == objsize2;
1138 /* iterate through key bits */
1139 for (i=0; i<depth; i++) {
1140 /* every second bit is a bit describing the object size */
1141 if (i%2 == 0) {
1142 /* check if object size bit is different in both keys (objsize1 and
1143 objsize2 describe the minimum index when object size bit is set) */
1144 if (
1145 (objsize1 <= j && objsize2 > j) ||
1146 (objsize2 <= j && objsize1 > j)
1147 ) {
1148 /* bit differs, therefore keys are in separate branches */
1149 return false;
1150 }
1151 /* increase bit counter for object size bits */
1152 j++;
1153 }
1154 /* all other bits describe latitude and longitude */
1155 else {
1156 /* check if bit differs in both keys */
1157 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
1158 /* bit differs, therefore keys are in separate branches */
1159 return false;
1160 }
1161 /* increase bit counter for latitude/longitude bits */
1162 k++;
1163 }
1164 }
1165 }
1166 /* if not, keys are point keys */
1167 else {
1168 /* iterate through key bits */
1169 for (i=0; i<depth; i++) {
1170 /* check if bit differs in both keys */
1171 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
1172 /* bit differs, therefore keys are in separate branches */
1173 return false;
1174 }
1175 }
1176 }
1177 /* return true because keys are in the same branch */
1178 return true;
1179 }
1181 /* combine two keys into new key which covers both original keys */
1182 /* (result stored in first argument) */
1183 static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
1184 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
1185 /* determine smallest depth */
1186 int depth1 = PGL_KEY_NODEDEPTH(dst);
1187 int depth2 = PGL_KEY_NODEDEPTH(src);
1188 int depth = (depth1 < depth2) ? depth1 : depth2;
1189 /* check if keys are area keys (assuming that both keys have same type) */
1190 if (PGL_KEY_IS_AREAKEY(dst)) {
1191 pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */
1192 int j = 0; /* bit offset for logarithmic object size bits */
1193 int k = 0; /* bit offset for latitude and longitude */
1194 /* fetch logarithmic object size information */
1195 int objsize1 = PGL_KEY_OBJSIZE(dst);
1196 int objsize2 = PGL_KEY_OBJSIZE(src);
1197 /* handle special cases for empty objects (universal and empty keys) */
1198 if (
1199 objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
1200 objsize2 > PGL_AREAKEY_MAXOBJSIZE
1201 ) {
1202 if (
1203 objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
1204 objsize2 == PGL_KEY_OBJSIZE_EMPTY
1205 ) pgl_key_set_empty(dst);
1206 else pgl_key_set_universal(dst);
1207 return;
1208 }
1209 /* iterate through key bits */
1210 for (i=0; i<depth; i++) {
1211 /* every second bit is a bit describing the object size */
1212 if (i%2 == 0) {
1213 /* increase bit counter for object size bits first */
1214 /* (handy when setting objsize variable) */
1215 j++;
1216 /* check if object size bit is set in neither key */
1217 if (objsize1 >= j && objsize2 >= j) {
1218 /* set objsize in destination buffer to indicate that size bit is
1219 unset in destination buffer at the current bit position */
1220 dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
1221 }
1222 /* break if object size bit is set in one key only */
1223 else if (objsize1 >= j || objsize2 >= j) break;
1224 }
1225 /* all other bits describe latitude and longitude */
1226 else {
1227 /* break if bit differs in both keys */
1228 if (PGL_KEY_LATLONBIT(dst, k)) {
1229 if (!PGL_KEY_LATLONBIT(src, k)) break;
1230 /* but set bit in destination buffer if bit is set in both keys */
1231 dstbuf[k/8] |= 0x80 >> (k%8);
1232 } else if (PGL_KEY_LATLONBIT(src, k)) break;
1233 /* increase bit counter for latitude/longitude bits */
1234 k++;
1235 }
1236 }
1237 /* set common node depth and type bit (type bit = 1) */
1238 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
1239 /* copy contents of destination buffer to first key */
1240 memcpy(dst, dstbuf, sizeof(pgl_areakey));
1241 }
1242 /* if not, keys are point keys */
1243 else {
1244 pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */
1245 /* iterate through key bits */
1246 for (i=0; i<depth; i++) {
1247 /* break if bit differs in both keys */
1248 if (PGL_KEY_LATLONBIT(dst, i)) {
1249 if (!PGL_KEY_LATLONBIT(src, i)) break;
1250 /* but set bit in destination buffer if bit is set in both keys */
1251 dstbuf[i/8] |= 0x80 >> (i%8);
1252 } else if (PGL_KEY_LATLONBIT(src, i)) break;
1253 }
1254 /* set common node depth (type bit = 0) */
1255 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
1256 /* copy contents of destination buffer to first key */
1257 memcpy(dst, dstbuf, sizeof(pgl_pointkey));
1258 }
1259 }
1261 /* determine center(!) boundaries and radius estimation of index key */
1262 static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
1263 int i;
1264 /* determine node depth */
1265 int depth = PGL_KEY_NODEDEPTH(key);
1266 /* center point of possible result */
1267 double lat = 0;
1268 double lon = 0;
1269 /* maximum distance of real center point from key center */
1270 double dlat = 90;
1271 double dlon = 180;
1272 /* maximum radius of contained objects */
1273 double radius = 0; /* always return zero for point index keys */
1274 /* check if key is area key */
1275 if (PGL_KEY_IS_AREAKEY(key)) {
1276 /* get logarithmic object size */
1277 int objsize = PGL_KEY_OBJSIZE(key);
1278 /* handle special cases for empty objects (universal and empty keys) */
1279 if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
1280 pgl_box_set_empty(box);
1281 return 0;
1282 } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
1283 box->lat_min = -90;
1284 box->lat_max = 90;
1285 box->lon_min = -180;
1286 box->lon_max = 180;
1287 return 0; /* any value >= 0 would do */
1288 }
1289 /* calculate maximum possible radius of objects covered by the given key */
1290 if (objsize == 0) radius = INFINITY;
1291 else {
1292 radius = PGL_AREAKEY_REFOBJSIZE;
1293 while (--objsize) radius /= M_SQRT2;
1294 }
1295 /* iterate over latitude and longitude bits in key */
1296 /* (every second bit is a latitude or longitude bit) */
1297 for (i=0; i<depth/2; i++) {
1298 /* check if latitude bit */
1299 if (i%2 == 0) {
1300 /* cut latitude dimension in half */
1301 dlat /= 2;
1302 /* increase center latitude if bit is 1, otherwise decrease */
1303 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1304 else lat -= dlat;
1305 }
1306 /* otherwise longitude bit */
1307 else {
1308 /* cut longitude dimension in half */
1309 dlon /= 2;
1310 /* increase center longitude if bit is 1, otherwise decrease */
1311 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1312 else lon -= dlon;
1313 }
1314 }
1315 }
1316 /* if not, keys are point keys */
1317 else {
1318 /* iterate over all bits in key */
1319 for (i=0; i<depth; i++) {
1320 /* check if latitude bit */
1321 if (i%2 == 0) {
1322 /* cut latitude dimension in half */
1323 dlat /= 2;
1324 /* increase center latitude if bit is 1, otherwise decrease */
1325 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1326 else lat -= dlat;
1327 }
1328 /* otherwise longitude bit */
1329 else {
1330 /* cut longitude dimension in half */
1331 dlon /= 2;
1332 /* increase center longitude if bit is 1, otherwise decrease */
1333 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1334 else lon -= dlon;
1335 }
1336 }
1337 }
1338 /* calculate boundaries from center point and remaining dlat and dlon */
1339 /* (return values through pointer to box) */
1340 box->lat_min = lat - dlat;
1341 box->lat_max = lat + dlat;
1342 box->lon_min = lon - dlon;
1343 box->lon_max = lon + dlon;
1344 /* return radius (as a function return value) */
1345 return radius;
1346 }
1348 /* estimator function for distance between point and index key */
1349 /* always returns a smaller value than actually correct or zero */
1350 static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
1351 pgl_box box; /* center(!) bounding box of area index key */
1352 /* calculate center(!) bounding box and maximum radius of objects covered
1353 by area index key (radius is zero for point index keys) */
1354 double distance = pgl_key_to_box(key, &box);
1355 /* calculate estimated distance between bounding box of center point of
1356 indexed object and point passed as second argument, then substract maximum
1357 radius of objects covered by index key */
1358 distance = pgl_estimate_point_box_distance(point, &box) - distance;
1359 /* truncate negative results to zero */
1360 if (distance <= 0) distance = 0;
1361 /* return result */
1362 return distance;
1363 }
1366 /*---------------------------------*
1367 * helper functions for text I/O *
1368 *---------------------------------*/
1370 #define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */
1372 /* convert floating point number to string (round-trip safe) */
1373 static void pgl_print_float(char *buf, double flt) {
1374 /* check if number is integral */
1375 if (trunc(flt) == flt) {
1376 /* for integral floats use maximum precision */
1377 snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1378 } else {
1379 /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
1380 snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
1381 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
1382 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1383 }
1384 }
1386 /* convert latitude floating point number (in degrees) to string */
1387 static void pgl_print_lat(char *buf, double lat) {
1388 if (signbit(lat)) {
1389 /* treat negative latitudes (including -0) as south */
1390 snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
1391 } else {
1392 /* treat positive latitudes (including +0) as north */
1393 snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
1394 }
1395 }
1397 /* convert longitude floating point number (in degrees) to string */
1398 static void pgl_print_lon(char *buf, double lon) {
1399 if (signbit(lon)) {
1400 /* treat negative longitudes (including -0) as west */
1401 snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
1402 } else {
1403 /* treat positive longitudes (including +0) as east */
1404 snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
1405 }
1406 }
1408 /* bit masks used as return value of pgl_scan() function */
1409 #define PGL_SCAN_NONE 0 /* no value has been parsed */
1410 #define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */
1411 #define PGL_SCAN_LON (1<<1) /* longitude has been parsed */
1412 #define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */
1414 /* parse a coordinate (can be latitude or longitude) */
1415 static int pgl_scan(char **str, double *lat, double *lon) {
1416 double val;
1417 int len;
1418 if (
1419 sscanf(*str, " N %lf %n", &val, &len) ||
1420 sscanf(*str, " n %lf %n", &val, &len)
1421 ) {
1422 *str += len; *lat = val; return PGL_SCAN_LAT;
1423 }
1424 if (
1425 sscanf(*str, " S %lf %n", &val, &len) ||
1426 sscanf(*str, " s %lf %n", &val, &len)
1427 ) {
1428 *str += len; *lat = -val; return PGL_SCAN_LAT;
1429 }
1430 if (
1431 sscanf(*str, " E %lf %n", &val, &len) ||
1432 sscanf(*str, " e %lf %n", &val, &len)
1433 ) {
1434 *str += len; *lon = val; return PGL_SCAN_LON;
1435 }
1436 if (
1437 sscanf(*str, " W %lf %n", &val, &len) ||
1438 sscanf(*str, " w %lf %n", &val, &len)
1439 ) {
1440 *str += len; *lon = -val; return PGL_SCAN_LON;
1441 }
1442 return PGL_SCAN_NONE;
1443 }
1446 /*-----------------*
1447 * SQL functions *
1448 *-----------------*/
1450 /* Note: These function names use "epoint", "ebox", etc. notation here instead
1451 of "point", "box", etc. in order to distinguish them from any previously
1452 defined functions. */
1454 /* function needed for dummy types and/or not implemented features */
1455 PG_FUNCTION_INFO_V1(pgl_notimpl);
1456 Datum pgl_notimpl(PG_FUNCTION_ARGS) {
1457 ereport(ERROR, (errmsg("not implemented by pgLatLon")));
1458 }
1460 /* set point to latitude and longitude (including checks) */
1461 static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
1462 /* reject infinite or NaN values */
1463 if (!isfinite(lat) || !isfinite(lon)) {
1464 ereport(ERROR, (
1465 errcode(ERRCODE_DATA_EXCEPTION),
1466 errmsg("epoint requires finite coordinates")
1467 ));
1468 }
1469 /* check latitude bounds */
1470 if (lat < -90) {
1471 ereport(WARNING, (errmsg("latitude exceeds south pole")));
1472 lat = -90;
1473 } else if (lat > 90) {
1474 ereport(WARNING, (errmsg("latitude exceeds north pole")));
1475 lat = 90;
1476 }
1477 /* check longitude bounds */
1478 if (lon < -180) {
1479 ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
1480 lon += 360 - trunc(lon / 360) * 360;
1481 } else if (lon > 180) {
1482 ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
1483 lon -= 360 + trunc(lon / 360) * 360;
1484 }
1485 /* store rounded latitude/longitude values for round-trip safety */
1486 point->lat = pgl_round(lat);
1487 point->lon = pgl_round(lon);
1488 }
1490 /* create point ("epoint" in SQL) from latitude and longitude */
1491 PG_FUNCTION_INFO_V1(pgl_create_epoint);
1492 Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
1493 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
1494 pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
1495 PG_RETURN_POINTER(point);
1496 }
1498 /* parse point ("epoint" in SQL) */
1499 /* format: '[NS]<float> [EW]<float>' */
1500 PG_FUNCTION_INFO_V1(pgl_epoint_in);
1501 Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
1502 char *str = PG_GETARG_CSTRING(0); /* input string */
1503 char *strptr = str; /* current position within string */
1504 int done = 0; /* bit mask storing if latitude or longitude was read */
1505 double lat, lon; /* parsed values as double precision floats */
1506 pgl_point *point; /* return value (to be palloc'ed) */
1507 /* parse two floats (each latitude or longitude) separated by white-space */
1508 done |= pgl_scan(&strptr, &lat, &lon);
1509 if (strptr != str && isspace(strptr[-1])) {
1510 done |= pgl_scan(&strptr, &lat, &lon);
1511 }
1512 /* require end of string, and latitude and longitude parsed successfully */
1513 if (strptr[0] || done != PGL_SCAN_LATLON) {
1514 ereport(ERROR, (
1515 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1516 errmsg("invalid input syntax for type epoint: \"%s\"", str)
1517 ));
1518 }
1519 /* allocate memory for result */
1520 point = (pgl_point *)palloc(sizeof(pgl_point));
1521 /* set latitude and longitude (and perform checks) */
1522 pgl_epoint_set_latlon(point, lat, lon);
1523 /* return result */
1524 PG_RETURN_POINTER(point);
1525 }
1527 /* create box ("ebox" in SQL) that is empty */
1528 PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
1529 Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
1530 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1531 pgl_box_set_empty(box);
1532 PG_RETURN_POINTER(box);
1533 }
1535 /* set box to given boundaries (including checks) */
1536 static void pgl_ebox_set_boundaries(
1537 pgl_box *box,
1538 double lat_min, double lat_max, double lon_min, double lon_max
1539 ) {
1540 /* if minimum latitude is greater than maximum latitude, return empty box */
1541 if (lat_min > lat_max) {
1542 pgl_box_set_empty(box);
1543 return;
1544 }
1545 /* otherwise reject infinite or NaN values */
1546 if (
1547 !isfinite(lat_min) || !isfinite(lat_max) ||
1548 !isfinite(lon_min) || !isfinite(lon_max)
1549 ) {
1550 ereport(ERROR, (
1551 errcode(ERRCODE_DATA_EXCEPTION),
1552 errmsg("ebox requires finite coordinates")
1553 ));
1554 }
1555 /* check latitude bounds */
1556 if (lat_max < -90) {
1557 ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
1558 lat_max = -90;
1559 } else if (lat_max > 90) {
1560 ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
1561 lat_max = 90;
1562 }
1563 if (lat_min < -90) {
1564 ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
1565 lat_min = -90;
1566 } else if (lat_min > 90) {
1567 ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
1568 lat_min = 90;
1569 }
1570 /* check if all longitudes are included */
1571 if (lon_max - lon_min >= 360) {
1572 if (lon_max - lon_min > 360) ereport(WARNING, (
1573 errmsg("longitude coverage greater than 360 degrees")
1574 ));
1575 lon_min = -180;
1576 lon_max = 180;
1577 } else {
1578 /* normalize longitude bounds */
1579 if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
1580 else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360;
1581 if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
1582 else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360;
1583 }
1584 /* store rounded latitude/longitude values for round-trip safety */
1585 box->lat_min = pgl_round(lat_min);
1586 box->lat_max = pgl_round(lat_max);
1587 box->lon_min = pgl_round(lon_min);
1588 box->lon_max = pgl_round(lon_max);
1589 /* ensure that rounding does not change orientation */
1590 if (lon_min > lon_max && box->lon_min == box->lon_max) {
1591 box->lon_min = -180;
1592 box->lon_max = 180;
1593 }
1594 }
1596 /* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
1597 PG_FUNCTION_INFO_V1(pgl_create_ebox);
1598 Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
1599 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1600 pgl_ebox_set_boundaries(
1601 box,
1602 PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
1603 PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
1604 );
1605 PG_RETURN_POINTER(box);
1606 }
1608 /* create box ("ebox" in SQL) from two points ("epoint"s) */
1609 /* (can not be used to cover a longitude range of more than 120 degrees) */
1610 PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
1611 Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
1612 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
1613 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
1614 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1615 double lat_min, lat_max, lon_min, lon_max;
1616 double dlon; /* longitude range (delta longitude) */
1617 /* order latitude and longitude boundaries */
1618 if (point2->lat < point1->lat) {
1619 lat_min = point2->lat;
1620 lat_max = point1->lat;
1621 } else {
1622 lat_min = point1->lat;
1623 lat_max = point2->lat;
1624 }
1625 if (point2->lon < point1->lon) {
1626 lon_min = point2->lon;
1627 lon_max = point1->lon;
1628 } else {
1629 lon_min = point1->lon;
1630 lon_max = point2->lon;
1631 }
1632 /* calculate longitude range (round to avoid floating point errors) */
1633 dlon = pgl_round(lon_max - lon_min);
1634 /* determine east-west direction */
1635 if (dlon >= 240) {
1636 /* assume that 180th meridian is crossed and swap min/max longitude */
1637 double swap = lon_min; lon_min = lon_max; lon_max = swap;
1638 } else if (dlon > 120) {
1639 /* unclear orientation since delta longitude > 120 */
1640 ereport(ERROR, (
1641 errcode(ERRCODE_DATA_EXCEPTION),
1642 errmsg("can not determine east/west orientation for ebox")
1643 ));
1644 }
1645 /* use boundaries to setup box (and perform checks) */
1646 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1647 /* return result */
1648 PG_RETURN_POINTER(box);
1649 }
1651 /* parse box ("ebox" in SQL) */
1652 /* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
1653 or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
1654 PG_FUNCTION_INFO_V1(pgl_ebox_in);
1655 Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
1656 char *str = PG_GETARG_CSTRING(0); /* input string */
1657 char *str_lower; /* lower case version of input string */
1658 char *strptr; /* current position within string */
1659 int valid; /* number of valid chars */
1660 int done; /* specifies if latitude or longitude was read */
1661 double val; /* temporary variable */
1662 int lat_count = 0; /* count of latitude values parsed */
1663 int lon_count = 0; /* count of longitufde values parsed */
1664 double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */
1665 pgl_box *box; /* return value (to be palloc'ed) */
1666 /* lowercase input */
1667 str_lower = psprintf("%s", str);
1668 for (strptr=str_lower; *strptr; strptr++) {
1669 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1670 }
1671 /* reset reading position to start of (lowercase) string */
1672 strptr = str_lower;
1673 /* check if empty box */
1674 valid = 0;
1675 sscanf(strptr, " empty %n", &valid);
1676 if (valid && strptr[valid] == 0) {
1677 /* allocate and return empty box */
1678 box = (pgl_box *)palloc(sizeof(pgl_box));
1679 pgl_box_set_empty(box);
1680 PG_RETURN_POINTER(box);
1681 }
1682 /* demand four blocks separated by whitespace */
1683 valid = 0;
1684 sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
1685 /* if four blocks separated by whitespace exist, parse those blocks */
1686 if (strptr[valid] == 0) while (strptr[0]) {
1687 /* parse either latitude or longitude (whichever found in input string) */
1688 done = pgl_scan(&strptr, &val, &val);
1689 /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
1690 if (done == PGL_SCAN_LAT) {
1691 if (!lat_count) lat_min = val; else lat_max = val;
1692 lat_count++;
1693 } else if (done == PGL_SCAN_LON) {
1694 if (!lon_count) lon_min = val; else lon_max = val;
1695 lon_count++;
1696 } else {
1697 break;
1698 }
1699 }
1700 /* require end of string, and two latitude and two longitude values */
1701 if (strptr[0] || lat_count != 2 || lon_count != 2) {
1702 ereport(ERROR, (
1703 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1704 errmsg("invalid input syntax for type ebox: \"%s\"", str)
1705 ));
1706 }
1707 /* free lower case string */
1708 pfree(str_lower);
1709 /* order boundaries (maximum greater than minimum) */
1710 if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
1711 if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
1712 /* allocate memory for result */
1713 box = (pgl_box *)palloc(sizeof(pgl_box));
1714 /* set boundaries (and perform checks) */
1715 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1716 /* return result */
1717 PG_RETURN_POINTER(box);
1718 }
1720 /* set circle to given latitude, longitude, and radius (including checks) */
1721 static void pgl_ecircle_set_latlon_radius(
1722 pgl_circle *circle, double lat, double lon, double radius
1723 ) {
1724 /* set center point (including checks) */
1725 pgl_epoint_set_latlon(&(circle->center), lat, lon);
1726 /* handle non-positive radius */
1727 if (isnan(radius)) {
1728 ereport(ERROR, (
1729 errcode(ERRCODE_DATA_EXCEPTION),
1730 errmsg("invalid radius for ecircle")
1731 ));
1732 }
1733 if (radius == 0) radius = 0; /* avoids -0 */
1734 else if (radius < 0) {
1735 if (isfinite(radius)) {
1736 ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
1737 }
1738 radius = -INFINITY;
1739 }
1740 /* store radius (round-trip safety is ensured by pgl_print_float) */
1741 circle->radius = radius;
1742 }
1744 /* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
1745 PG_FUNCTION_INFO_V1(pgl_create_ecircle);
1746 Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
1747 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1748 pgl_ecircle_set_latlon_radius(
1749 circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
1750 );
1751 PG_RETURN_POINTER(circle);
1752 }
1754 /* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
1755 PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
1756 Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
1757 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1758 double radius = PG_GETARG_FLOAT8(1);
1759 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1760 /* set latitude, longitude, radius (and perform checks) */
1761 pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
1762 /* return result */
1763 PG_RETURN_POINTER(circle);
1764 }
1766 /* parse circle ("ecircle" in SQL) */
1767 /* format: '[NS]<float> [EW]<float> <float>' */
1768 PG_FUNCTION_INFO_V1(pgl_ecircle_in);
1769 Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
1770 char *str = PG_GETARG_CSTRING(0); /* input string */
1771 char *strptr = str; /* current position within string */
1772 double lat, lon, radius; /* parsed values as double precision flaots */
1773 int valid = 0; /* number of valid chars */
1774 int done = 0; /* stores if latitude and/or longitude was read */
1775 pgl_circle *circle; /* return value (to be palloc'ed) */
1776 /* demand three blocks separated by whitespace */
1777 sscanf(strptr, " %*s %*s %*s %n", &valid);
1778 /* if three blocks separated by whitespace exist, parse those blocks */
1779 if (strptr[valid] == 0) {
1780 /* parse latitude and longitude */
1781 done |= pgl_scan(&strptr, &lat, &lon);
1782 done |= pgl_scan(&strptr, &lat, &lon);
1783 /* parse radius (while incrementing strptr by number of bytes parsed) */
1784 valid = 0;
1785 if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
1786 }
1787 /* require end of string and both latitude and longitude being parsed */
1788 if (strptr[0] || done != PGL_SCAN_LATLON) {
1789 ereport(ERROR, (
1790 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1791 errmsg("invalid input syntax for type ecircle: \"%s\"", str)
1792 ));
1793 }
1794 /* allocate memory for result */
1795 circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1796 /* set latitude, longitude, radius (and perform checks) */
1797 pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
1798 /* return result */
1799 PG_RETURN_POINTER(circle);
1800 }
1802 /* parse cluster ("ecluster" in SQL) */
1803 PG_FUNCTION_INFO_V1(pgl_ecluster_in);
1804 Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
1805 int i;
1806 char *str = PG_GETARG_CSTRING(0); /* input string */
1807 char *str_lower; /* lower case version of input string */
1808 char *strptr; /* pointer to current reading position of input */
1809 int npoints_total = 0; /* total number of points in cluster */
1810 int nentries = 0; /* total number of entries */
1811 pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */
1812 int entries_buflen = 4; /* maximum number of elements in entries array */
1813 int valid; /* number of valid chars processed */
1814 double lat, lon; /* latitude and longitude of parsed point */
1815 int entrytype; /* current entry type */
1816 int npoints; /* number of points in current entry */
1817 pgl_point *points; /* array of pgl_point for pgl_newentry */
1818 int points_buflen; /* maximum number of elements in points array */
1819 int done; /* return value of pgl_scan function */
1820 pgl_cluster *cluster; /* created cluster */
1821 /* lowercase input */
1822 str_lower = psprintf("%s", str);
1823 for (strptr=str_lower; *strptr; strptr++) {
1824 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1825 }
1826 /* reset reading position to start of (lowercase) string */
1827 strptr = str_lower;
1828 /* allocate initial buffer for entries */
1829 entries = palloc(entries_buflen * sizeof(pgl_newentry));
1830 /* parse until end of string */
1831 while (strptr[0]) {
1832 /* require previous white-space or closing parenthesis before next token */
1833 if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
1834 goto pgl_ecluster_in_error;
1835 }
1836 /* ignore token "empty" */
1837 valid = 0; sscanf(strptr, " empty %n", &valid);
1838 if (valid) { strptr += valid; continue; }
1839 /* test for "point" token */
1840 valid = 0; sscanf(strptr, " point ( %n", &valid);
1841 if (valid) {
1842 strptr += valid;
1843 entrytype = PGL_ENTRY_POINT;
1844 goto pgl_ecluster_in_type_ok;
1845 }
1846 /* test for "path" token */
1847 valid = 0; sscanf(strptr, " path ( %n", &valid);
1848 if (valid) {
1849 strptr += valid;
1850 entrytype = PGL_ENTRY_PATH;
1851 goto pgl_ecluster_in_type_ok;
1852 }
1853 /* test for "outline" token */
1854 valid = 0; sscanf(strptr, " outline ( %n", &valid);
1855 if (valid) {
1856 strptr += valid;
1857 entrytype = PGL_ENTRY_OUTLINE;
1858 goto pgl_ecluster_in_type_ok;
1859 }
1860 /* test for "polygon" token */
1861 valid = 0; sscanf(strptr, " polygon ( %n", &valid);
1862 if (valid) {
1863 strptr += valid;
1864 entrytype = PGL_ENTRY_POLYGON;
1865 goto pgl_ecluster_in_type_ok;
1866 }
1867 /* error if no valid token found */
1868 goto pgl_ecluster_in_error;
1869 pgl_ecluster_in_type_ok:
1870 /* check if pgl_newentry array needs to grow */
1871 if (nentries == entries_buflen) {
1872 pgl_newentry *newbuf;
1873 entries_buflen *= 2;
1874 newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
1875 memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
1876 pfree(entries);
1877 entries = newbuf;
1878 }
1879 /* reset number of points for current entry */
1880 npoints = 0;
1881 /* allocate array for points */
1882 points_buflen = 4;
1883 points = palloc(points_buflen * sizeof(pgl_point));
1884 /* parse until closing parenthesis */
1885 while (strptr[0] != ')') {
1886 /* error on unexpected end of string */
1887 if (strptr[0] == 0) goto pgl_ecluster_in_error;
1888 /* mark neither latitude nor longitude as read */
1889 done = PGL_SCAN_NONE;
1890 /* require white-space before second, third, etc. point */
1891 if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1892 /* scan latitude (or longitude) */
1893 done |= pgl_scan(&strptr, &lat, &lon);
1894 /* require white-space before second coordinate */
1895 if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1896 /* scan longitude (or latitude) */
1897 done |= pgl_scan(&strptr, &lat, &lon);
1898 /* error unless both latitude and longitude were parsed */
1899 if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
1900 /* throw error if number of points is too high */
1901 if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
1902 ereport(ERROR, (
1903 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1904 errmsg(
1905 "too many points for ecluster entry (maximum %i)",
1906 PGL_CLUSTER_MAXPOINTS
1907 )
1908 ));
1909 }
1910 /* check if pgl_point array needs to grow */
1911 if (npoints == points_buflen) {
1912 pgl_point *newbuf;
1913 points_buflen *= 2;
1914 newbuf = palloc(points_buflen * sizeof(pgl_point));
1915 memcpy(newbuf, points, npoints * sizeof(pgl_point));
1916 pfree(points);
1917 points = newbuf;
1918 }
1919 /* append point to pgl_point array (includes checks) */
1920 pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
1921 /* increase total number of points */
1922 npoints_total++;
1923 }
1924 /* error if entry has no points */
1925 if (!npoints) goto pgl_ecluster_in_error;
1926 /* entries with one point are automatically of type "point" */
1927 if (npoints == 1) entrytype = PGL_ENTRY_POINT;
1928 /* if entries have more than one point */
1929 else {
1930 /* throw error if entry type is "point" */
1931 if (entrytype == PGL_ENTRY_POINT) {
1932 ereport(ERROR, (
1933 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1934 errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
1935 ));
1936 }
1937 /* coerce outlines and polygons with more than 2 points to be a path */
1938 if (npoints == 2) entrytype = PGL_ENTRY_PATH;
1939 }
1940 /* append entry to pgl_newentry array */
1941 entries[nentries].entrytype = entrytype;
1942 entries[nentries].npoints = npoints;
1943 entries[nentries].points = points;
1944 nentries++;
1945 /* consume closing parenthesis */
1946 strptr++;
1947 /* consume white-space */
1948 while (isspace(strptr[0])) strptr++;
1949 }
1950 /* free lower case string */
1951 pfree(str_lower);
1952 /* create cluster from pgl_newentry array */
1953 cluster = pgl_new_cluster(nentries, entries);
1954 /* free pgl_newentry array */
1955 for (i=0; i<nentries; i++) pfree(entries[i].points);
1956 pfree(entries);
1957 /* set bounding circle of cluster and check east/west orientation */
1958 if (!pgl_finalize_cluster(cluster)) {
1959 ereport(ERROR, (
1960 errcode(ERRCODE_DATA_EXCEPTION),
1961 errmsg("can not determine east/west orientation for ecluster"),
1962 errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
1963 ));
1964 }
1965 /* return cluster */
1966 PG_RETURN_POINTER(cluster);
1967 /* code to throw error */
1968 pgl_ecluster_in_error:
1969 ereport(ERROR, (
1970 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1971 errmsg("invalid input syntax for type ecluster: \"%s\"", str)
1972 ));
1973 }
1975 /* convert point ("epoint") to string representation */
1976 PG_FUNCTION_INFO_V1(pgl_epoint_out);
1977 Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
1978 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1979 char latstr[PGL_NUMBUFLEN];
1980 char lonstr[PGL_NUMBUFLEN];
1981 pgl_print_lat(latstr, point->lat);
1982 pgl_print_lon(lonstr, point->lon);
1983 PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
1984 }
1986 /* convert box ("ebox") to string representation */
1987 PG_FUNCTION_INFO_V1(pgl_ebox_out);
1988 Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
1989 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
1990 double lon_min = box->lon_min;
1991 double lon_max = box->lon_max;
1992 char lat_min_str[PGL_NUMBUFLEN];
1993 char lat_max_str[PGL_NUMBUFLEN];
1994 char lon_min_str[PGL_NUMBUFLEN];
1995 char lon_max_str[PGL_NUMBUFLEN];
1996 /* return string "empty" if box is set to be empty */
1997 if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
1998 /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
1999 /* (required since pgl_box_in orders the longitude boundaries) */
2000 if (lon_min > lon_max) {
2001 if (lon_min + lon_max >= 0) lon_min -= 360;
2002 else lon_max += 360;
2003 }
2004 /* format and return result */
2005 pgl_print_lat(lat_min_str, box->lat_min);
2006 pgl_print_lat(lat_max_str, box->lat_max);
2007 pgl_print_lon(lon_min_str, lon_min);
2008 pgl_print_lon(lon_max_str, lon_max);
2009 PG_RETURN_CSTRING(psprintf(
2010 "%s %s %s %s",
2011 lat_min_str, lon_min_str, lat_max_str, lon_max_str
2012 ));
2013 }
2015 /* convert circle ("ecircle") to string representation */
2016 PG_FUNCTION_INFO_V1(pgl_ecircle_out);
2017 Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
2018 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2019 char latstr[PGL_NUMBUFLEN];
2020 char lonstr[PGL_NUMBUFLEN];
2021 char radstr[PGL_NUMBUFLEN];
2022 pgl_print_lat(latstr, circle->center.lat);
2023 pgl_print_lon(lonstr, circle->center.lon);
2024 pgl_print_float(radstr, circle->radius);
2025 PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
2026 }
2028 /* convert cluster ("ecluster") to string representation */
2029 PG_FUNCTION_INFO_V1(pgl_ecluster_out);
2030 Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
2031 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2032 char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */
2033 char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */
2034 char ***strings; /* array of array of strings */
2035 char *string; /* string of current token */
2036 char *res, *resptr; /* result and pointer to current write position */
2037 size_t reslen = 1; /* length of result (init with 1 for terminator) */
2038 int npoints; /* number of points of current entry */
2039 int i, j; /* i: entry, j: point in entry */
2040 /* handle empty clusters */
2041 if (cluster->nentries == 0) {
2042 /* free detoasted cluster (if copy) */
2043 PG_FREE_IF_COPY(cluster, 0);
2044 /* return static result */
2045 PG_RETURN_CSTRING("empty");
2046 }
2047 /* allocate array of array of strings */
2048 strings = palloc(cluster->nentries * sizeof(char **));
2049 /* iterate over all entries in cluster */
2050 for (i=0; i<cluster->nentries; i++) {
2051 /* get number of points in entry */
2052 npoints = cluster->entries[i].npoints;
2053 /* allocate array of strings (one string for each point plus two extra) */
2054 strings[i] = palloc((2 + npoints) * sizeof(char *));
2055 /* determine opening string */
2056 switch (cluster->entries[i].entrytype) {
2057 case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break;
2058 case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break;
2059 case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
2060 case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
2061 default: string = (i==0)?"unknown" :" unknown";
2062 }
2063 /* use opening string as first string in array */
2064 strings[i][0] = string;
2065 /* update result length (for allocating result string later) */
2066 reslen += strlen(string);
2067 /* iterate over all points */
2068 for (j=0; j<npoints; j++) {
2069 /* create string representation of point */
2070 pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
2071 pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
2072 string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
2073 /* copy string pointer to string array */
2074 strings[i][j+1] = string;
2075 /* update result length (for allocating result string later) */
2076 reslen += strlen(string);
2077 }
2078 /* use closing parenthesis as last string in array */
2079 strings[i][npoints+1] = ")";
2080 /* update result length (for allocating result string later) */
2081 reslen++;
2082 }
2083 /* allocate result string */
2084 res = palloc(reslen);
2085 /* set write pointer to begin of result string */
2086 resptr = res;
2087 /* copy strings into result string */
2088 for (i=0; i<cluster->nentries; i++) {
2089 npoints = cluster->entries[i].npoints;
2090 for (j=0; j<npoints+2; j++) {
2091 string = strings[i][j];
2092 strcpy(resptr, string);
2093 resptr += strlen(string);
2094 /* free strings allocated by psprintf */
2095 if (j != 0 && j != npoints+1) pfree(string);
2096 }
2097 /* free array of strings */
2098 pfree(strings[i]);
2099 }
2100 /* free array of array of strings */
2101 pfree(strings);
2102 /* free detoasted cluster (if copy) */
2103 PG_FREE_IF_COPY(cluster, 0);
2104 /* return result */
2105 PG_RETURN_CSTRING(res);
2106 }
2108 /* binary input function for point ("epoint") */
2109 PG_FUNCTION_INFO_V1(pgl_epoint_recv);
2110 Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
2111 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2112 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
2113 point->lat = pq_getmsgfloat8(buf);
2114 point->lon = pq_getmsgfloat8(buf);
2115 PG_RETURN_POINTER(point);
2116 }
2118 /* binary input function for box ("ebox") */
2119 PG_FUNCTION_INFO_V1(pgl_ebox_recv);
2120 Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
2121 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2122 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
2123 box->lat_min = pq_getmsgfloat8(buf);
2124 box->lat_max = pq_getmsgfloat8(buf);
2125 box->lon_min = pq_getmsgfloat8(buf);
2126 box->lon_max = pq_getmsgfloat8(buf);
2127 PG_RETURN_POINTER(box);
2128 }
2130 /* binary input function for circle ("ecircle") */
2131 PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
2132 Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
2133 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
2134 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
2135 circle->center.lat = pq_getmsgfloat8(buf);
2136 circle->center.lon = pq_getmsgfloat8(buf);
2137 circle->radius = pq_getmsgfloat8(buf);
2138 PG_RETURN_POINTER(circle);
2139 }
2141 /* TODO: binary receive function for cluster */
2143 /* binary output function for point ("epoint") */
2144 PG_FUNCTION_INFO_V1(pgl_epoint_send);
2145 Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
2146 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2147 StringInfoData buf;
2148 pq_begintypsend(&buf);
2149 pq_sendfloat8(&buf, point->lat);
2150 pq_sendfloat8(&buf, point->lon);
2151 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2152 }
2154 /* binary output function for box ("ebox") */
2155 PG_FUNCTION_INFO_V1(pgl_ebox_send);
2156 Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
2157 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2158 StringInfoData buf;
2159 pq_begintypsend(&buf);
2160 pq_sendfloat8(&buf, box->lat_min);
2161 pq_sendfloat8(&buf, box->lat_max);
2162 pq_sendfloat8(&buf, box->lon_min);
2163 pq_sendfloat8(&buf, box->lon_max);
2164 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2165 }
2167 /* binary output function for circle ("ecircle") */
2168 PG_FUNCTION_INFO_V1(pgl_ecircle_send);
2169 Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
2170 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2171 StringInfoData buf;
2172 pq_begintypsend(&buf);
2173 pq_sendfloat8(&buf, circle->center.lat);
2174 pq_sendfloat8(&buf, circle->center.lon);
2175 pq_sendfloat8(&buf, circle->radius);
2176 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2177 }
2179 /* TODO: binary send functions for cluster */
2181 /* cast point ("epoint") to box ("ebox") */
2182 PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
2183 Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
2184 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2185 pgl_box *box = palloc(sizeof(pgl_box));
2186 box->lat_min = point->lat;
2187 box->lat_max = point->lat;
2188 box->lon_min = point->lon;
2189 box->lon_max = point->lon;
2190 PG_RETURN_POINTER(box);
2191 }
2193 /* cast point ("epoint") to circle ("ecircle") */
2194 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
2195 Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
2196 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2197 pgl_circle *circle = palloc(sizeof(pgl_box));
2198 circle->center = *point;
2199 circle->radius = 0;
2200 PG_RETURN_POINTER(circle);
2201 }
2203 /* cast point ("epoint") to cluster ("ecluster") */
2204 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
2205 Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
2206 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2207 pgl_newentry entry;
2208 entry.entrytype = PGL_ENTRY_POINT;
2209 entry.npoints = 1;
2210 entry.points = point;
2211 PG_RETURN_POINTER(pgl_new_cluster(1, &entry));
2212 }
2214 /* cast box ("ebox") to cluster ("ecluster") */
2215 #define pgl_ebox_to_ecluster_macro(i, a, b) \
2216 entries[i].entrytype = PGL_ENTRY_POLYGON; \
2217 entries[i].npoints = 4; \
2218 entries[i].points = points[i]; \
2219 points[i][0].lat = box->lat_min; \
2220 points[i][0].lon = (a); \
2221 points[i][1].lat = box->lat_min; \
2222 points[i][1].lon = (b); \
2223 points[i][2].lat = box->lat_max; \
2224 points[i][2].lon = (b); \
2225 points[i][3].lat = box->lat_max; \
2226 points[i][3].lon = (a);
2227 PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
2228 Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
2229 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2230 double lon, dlon;
2231 int nentries;
2232 pgl_newentry entries[3];
2233 pgl_point points[3][4];
2234 if (box->lat_min > box->lat_max) {
2235 nentries = 0;
2236 } else if (box->lon_min > box->lon_max) {
2237 if (box->lon_min < 0) {
2238 lon = pgl_round((box->lon_min + 180) / 2.0);
2239 nentries = 3;
2240 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2241 pgl_ebox_to_ecluster_macro(1, lon, 180);
2242 pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
2243 } else if (box->lon_max > 0) {
2244 lon = pgl_round((box->lon_max - 180) / 2.0);
2245 nentries = 3;
2246 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2247 pgl_ebox_to_ecluster_macro(1, -180, lon);
2248 pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
2249 } else {
2250 nentries = 2;
2251 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2252 pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
2253 }
2254 } else {
2255 dlon = pgl_round(box->lon_max - box->lon_min);
2256 if (dlon < 180) {
2257 nentries = 1;
2258 pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
2259 } else {
2260 lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
2261 if (
2262 pgl_round(lon - box->lon_min) < 180 &&
2263 pgl_round(box->lon_max - lon) < 180
2264 ) {
2265 nentries = 2;
2266 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2267 pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
2268 } else {
2269 nentries = 3;
2270 pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
2271 pgl_ebox_to_ecluster_macro(1, -60, 60);
2272 pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
2273 }
2274 }
2275 }
2276 PG_RETURN_POINTER(pgl_new_cluster(nentries, entries));
2277 }
2279 /* extract latitude from point ("epoint") */
2280 PG_FUNCTION_INFO_V1(pgl_epoint_lat);
2281 Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
2282 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
2283 }
2285 /* extract longitude from point ("epoint") */
2286 PG_FUNCTION_INFO_V1(pgl_epoint_lon);
2287 Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
2288 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
2289 }
2291 /* extract minimum latitude from box ("ebox") */
2292 PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
2293 Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
2294 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
2295 }
2297 /* extract maximum latitude from box ("ebox") */
2298 PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
2299 Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
2300 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
2301 }
2303 /* extract minimum longitude from box ("ebox") */
2304 PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
2305 Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
2306 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
2307 }
2309 /* extract maximum longitude from box ("ebox") */
2310 PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
2311 Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
2312 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
2313 }
2315 /* extract center point from circle ("ecircle") */
2316 PG_FUNCTION_INFO_V1(pgl_ecircle_center);
2317 Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
2318 PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
2319 }
2321 /* extract radius from circle ("ecircle") */
2322 PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
2323 Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
2324 PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
2325 }
2327 /* check if point is inside box (overlap operator "&&") in SQL */
2328 PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
2329 Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
2330 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2331 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
2332 PG_RETURN_BOOL(pgl_point_in_box(point, box));
2333 }
2335 /* check if point is inside circle (overlap operator "&&") in SQL */
2336 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
2337 Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
2338 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2339 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2340 PG_RETURN_BOOL(
2341 pgl_distance(
2342 point->lat, point->lon,
2343 circle->center.lat, circle->center.lon
2344 ) <= circle->radius
2345 );
2346 }
2348 /* check if point is inside cluster (overlap operator "&&") in SQL */
2349 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
2350 Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
2351 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2352 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2353 bool retval;
2354 /* points outside bounding circle are always assumed to be non-overlapping
2355 (necessary for consistent table and index scans) */
2356 if (
2357 pgl_distance(
2358 point->lat, point->lon,
2359 cluster->bounding.center.lat, cluster->bounding.center.lon
2360 ) > cluster->bounding.radius
2361 ) retval = false;
2362 else retval = pgl_point_in_cluster(point, cluster, false);
2363 PG_FREE_IF_COPY(cluster, 1);
2364 PG_RETURN_BOOL(retval);
2365 }
2367 /* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
2368 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
2369 Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2370 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2371 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2372 bool retval = pgl_distance(
2373 point->lat, point->lon,
2374 cluster->bounding.center.lat, cluster->bounding.center.lon
2375 ) <= cluster->bounding.radius;
2376 PG_FREE_IF_COPY(cluster, 1);
2377 PG_RETURN_BOOL(retval);
2378 }
2380 /* check if two boxes overlap (overlap operator "&&") in SQL */
2381 PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
2382 Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
2383 pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
2384 pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
2385 PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
2386 }
2388 /* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
2389 PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
2390 Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
2391 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2392 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2393 PG_RETURN_BOOL(
2394 pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
2395 );
2396 }
2398 /* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
2399 PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
2400 Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2401 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
2402 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2403 bool retval = pgl_estimate_point_box_distance(
2404 &cluster->bounding.center,
2405 box
2406 ) <= cluster->bounding.radius;
2407 PG_FREE_IF_COPY(cluster, 1);
2408 PG_RETURN_BOOL(retval);
2409 }
2411 /* check if two circles overlap (overlap operator "&&") in SQL */
2412 PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
2413 Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
2414 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2415 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2416 PG_RETURN_BOOL(
2417 pgl_distance(
2418 circle1->center.lat, circle1->center.lon,
2419 circle2->center.lat, circle2->center.lon
2420 ) <= circle1->radius + circle2->radius
2421 );
2422 }
2424 /* check if circle and cluster overlap (overlap operator "&&") in SQL */
2425 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
2426 Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
2427 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2428 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2429 bool retval = (
2430 pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
2431 );
2432 PG_FREE_IF_COPY(cluster, 1);
2433 PG_RETURN_BOOL(retval);
2434 }
2436 /* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
2437 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
2438 Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2439 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2440 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2441 bool retval = pgl_distance(
2442 circle->center.lat, circle->center.lon,
2443 cluster->bounding.center.lat, cluster->bounding.center.lon
2444 ) <= circle->radius + cluster->bounding.radius;
2445 PG_FREE_IF_COPY(cluster, 1);
2446 PG_RETURN_BOOL(retval);
2447 }
2449 /* check if two clusters overlap (overlap operator "&&") in SQL */
2450 PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
2451 Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
2452 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2453 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2454 bool retval;
2455 /* clusters with non-touching bounding circles are always assumed to be
2456 non-overlapping (improves performance and is necessary for consistent
2457 table and index scans) */
2458 if (
2459 pgl_distance(
2460 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
2461 cluster2->bounding.center.lat, cluster2->bounding.center.lon
2462 ) > cluster1->bounding.radius + cluster2->bounding.radius
2463 ) retval = false;
2464 else retval = pgl_clusters_overlap(cluster1, cluster2);
2465 PG_FREE_IF_COPY(cluster1, 0);
2466 PG_FREE_IF_COPY(cluster2, 1);
2467 PG_RETURN_BOOL(retval);
2468 }
2470 /* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
2471 PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
2472 Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
2473 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2474 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2475 bool retval = pgl_distance(
2476 cluster1->bounding.center.lat, cluster1->bounding.center.lon,
2477 cluster2->bounding.center.lat, cluster2->bounding.center.lon
2478 ) <= cluster1->bounding.radius + cluster2->bounding.radius;
2479 PG_FREE_IF_COPY(cluster1, 0);
2480 PG_FREE_IF_COPY(cluster2, 1);
2481 PG_RETURN_BOOL(retval);
2482 }
2484 /* check if second cluster is in first cluster (cont. operator "@>) in SQL */
2485 PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
2486 Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
2487 pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2488 pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2489 bool retval;
2490 /* clusters with non-touching bounding circles are always assumed to be
2491 non-overlapping (improves performance and is necessary for consistent
2492 table and index scans) */
2493 if (
2494 pgl_distance(
2495 outer->bounding.center.lat, outer->bounding.center.lon,
2496 inner->bounding.center.lat, inner->bounding.center.lon
2497 ) > outer->bounding.radius + inner->bounding.radius
2498 ) retval = false;
2499 else retval = pgl_cluster_in_cluster(outer, inner);
2500 PG_FREE_IF_COPY(outer, 0);
2501 PG_FREE_IF_COPY(inner, 1);
2502 PG_RETURN_BOOL(retval);
2503 }
2505 /* calculate distance between two points ("<->" operator) in SQL */
2506 PG_FUNCTION_INFO_V1(pgl_epoint_distance);
2507 Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
2508 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
2509 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
2510 PG_RETURN_FLOAT8(pgl_distance(
2511 point1->lat, point1->lon, point2->lat, point2->lon
2512 ));
2513 }
2515 /* calculate point to circle distance ("<->" operator) in SQL */
2516 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
2517 Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
2518 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2519 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2520 double distance = pgl_distance(
2521 point->lat, point->lon, circle->center.lat, circle->center.lon
2522 ) - circle->radius;
2523 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2524 }
2526 /* calculate point to cluster distance ("<->" operator) in SQL */
2527 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
2528 Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
2529 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2530 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2531 double distance = pgl_point_cluster_distance(point, cluster);
2532 PG_FREE_IF_COPY(cluster, 1);
2533 PG_RETURN_FLOAT8(distance);
2534 }
2536 /* calculate distance between two circles ("<->" operator) in SQL */
2537 PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
2538 Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
2539 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2540 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2541 double distance = pgl_distance(
2542 circle1->center.lat, circle1->center.lon,
2543 circle2->center.lat, circle2->center.lon
2544 ) - (circle1->radius + circle2->radius);
2545 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2546 }
2548 /* calculate circle to cluster distance ("<->" operator) in SQL */
2549 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
2550 Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
2551 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2552 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2553 double distance = (
2554 pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
2555 );
2556 PG_FREE_IF_COPY(cluster, 1);
2557 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2558 }
2560 /* calculate distance between two clusters ("<->" operator) in SQL */
2561 PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
2562 Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
2563 pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2564 pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2565 double retval = pgl_cluster_distance(cluster1, cluster2);
2566 PG_FREE_IF_COPY(cluster1, 0);
2567 PG_FREE_IF_COPY(cluster2, 1);
2568 PG_RETURN_FLOAT8(retval);
2569 }
2572 /*-----------------------------------------------------------*
2573 * B-tree comparison operators and index support functions *
2574 *-----------------------------------------------------------*/
2576 /* macro for a B-tree operator (without detoasting) */
2577 #define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
2578 PG_FUNCTION_INFO_V1(func); \
2579 Datum func(PG_FUNCTION_ARGS) { \
2580 type *a = (type *)PG_GETARG_POINTER(0); \
2581 type *b = (type *)PG_GETARG_POINTER(1); \
2582 PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
2583 }
2585 /* macro for a B-tree comparison function (without detoasting) */
2586 #define PGL_BTREE_CMP(func, type, cmpfunc) \
2587 PG_FUNCTION_INFO_V1(func); \
2588 Datum func(PG_FUNCTION_ARGS) { \
2589 type *a = (type *)PG_GETARG_POINTER(0); \
2590 type *b = (type *)PG_GETARG_POINTER(1); \
2591 PG_RETURN_INT32(cmpfunc(a, b)); \
2592 }
2594 /* macro for a B-tree operator (with detoasting) */
2595 #define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
2596 PG_FUNCTION_INFO_V1(func); \
2597 Datum func(PG_FUNCTION_ARGS) { \
2598 bool res; \
2599 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2600 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2601 res = cmpfunc(a, b) oper 0; \
2602 PG_FREE_IF_COPY(a, 0); \
2603 PG_FREE_IF_COPY(b, 1); \
2604 PG_RETURN_BOOL(res); \
2605 }
2607 /* macro for a B-tree comparison function (with detoasting) */
2608 #define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
2609 PG_FUNCTION_INFO_V1(func); \
2610 Datum func(PG_FUNCTION_ARGS) { \
2611 int32_t res; \
2612 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2613 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2614 res = cmpfunc(a, b); \
2615 PG_FREE_IF_COPY(a, 0); \
2616 PG_FREE_IF_COPY(b, 1); \
2617 PG_RETURN_INT32(res); \
2618 }
2620 /* B-tree operators and comparison function for point */
2621 PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
2622 PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
2623 PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
2624 PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
2625 PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
2626 PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
2627 PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
2629 /* B-tree operators and comparison function for box */
2630 PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
2631 PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
2632 PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
2633 PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
2634 PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
2635 PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
2636 PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
2638 /* B-tree operators and comparison function for circle */
2639 PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
2640 PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
2641 PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
2642 PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
2643 PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
2644 PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
2645 PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
2648 /*--------------------------------*
2649 * GiST index support functions *
2650 *--------------------------------*/
2652 /* GiST "consistent" support function */
2653 PG_FUNCTION_INFO_V1(pgl_gist_consistent);
2654 Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
2655 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2656 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
2657 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
2658 bool *recheck = (bool *)PG_GETARG_POINTER(4);
2659 /* demand recheck because index and query methods are lossy */
2660 *recheck = true;
2661 /* strategy number aliases for different operators using the same strategy */
2662 strategy %= 100;
2663 /* strategy number 11: equality of two points */
2664 if (strategy == 11) {
2665 /* query datum is another point */
2666 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2667 /* convert other point to key */
2668 pgl_pointkey querykey;
2669 pgl_point_to_key(query, querykey);
2670 /* return true if both keys overlap */
2671 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2672 }
2673 /* strategy number 13: equality of two circles */
2674 if (strategy == 13) {
2675 /* query datum is another circle */
2676 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2677 /* convert other circle to key */
2678 pgl_areakey querykey;
2679 pgl_circle_to_key(query, querykey);
2680 /* return true if both keys overlap */
2681 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2682 }
2683 /* for all remaining strategies, keys on empty objects produce no match */
2684 /* (check necessary because query radius may be infinite) */
2685 if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
2686 /* strategy number 21: overlapping with point */
2687 if (strategy == 21) {
2688 /* query datum is a point */
2689 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2690 /* return true if estimated distance (allowed to be smaller than real
2691 distance) between index key and point is zero */
2692 PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
2693 }
2694 /* strategy number 22: (point) overlapping with box */
2695 if (strategy == 22) {
2696 /* query datum is a box */
2697 pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
2698 /* determine bounding box of indexed key */
2699 pgl_box keybox;
2700 pgl_key_to_box(key, &keybox);
2701 /* return true if query box overlaps with bounding box of indexed key */
2702 PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
2703 }
2704 /* strategy number 23: overlapping with circle */
2705 if (strategy == 23) {
2706 /* query datum is a circle */
2707 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2708 /* return true if estimated distance (allowed to be smaller than real
2709 distance) between index key and circle center is smaller than radius */
2710 PG_RETURN_BOOL(
2711 pgl_estimate_key_distance(key, &(query->center)) <= query->radius
2712 );
2713 }
2714 /* strategy number 24: overlapping with cluster */
2715 if (strategy == 24) {
2716 bool retval; /* return value */
2717 /* query datum is a cluster */
2718 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2719 /* return true if estimated distance (allowed to be smaller than real
2720 distance) between index key and circle center is smaller than radius */
2721 retval = (
2722 pgl_estimate_key_distance(key, &(query->bounding.center)) <=
2723 query->bounding.radius
2724 );
2725 PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */
2726 PG_RETURN_BOOL(retval);
2727 }
2728 /* throw error for any unknown strategy number */
2729 elog(ERROR, "unrecognized strategy number: %d", strategy);
2730 }
2732 /* GiST "union" support function */
2733 PG_FUNCTION_INFO_V1(pgl_gist_union);
2734 Datum pgl_gist_union(PG_FUNCTION_ARGS) {
2735 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2736 pgl_keyptr out; /* return value (to be palloc'ed) */
2737 int i;
2738 /* determine key size */
2739 size_t keysize = PGL_KEY_IS_AREAKEY(
2740 (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
2741 ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
2742 /* begin with first key as result */
2743 out = palloc(keysize);
2744 memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
2745 /* unite current result with second, third, etc. key */
2746 for (i=1; i<entryvec->n; i++) {
2747 pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
2748 }
2749 /* return result */
2750 PG_RETURN_POINTER(out);
2751 }
2753 /* GiST "compress" support function for indicis on points */
2754 PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
2755 Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
2756 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2757 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2758 /* only transform new leaves */
2759 if (entry->leafkey) {
2760 /* get point to be transformed */
2761 pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
2762 /* allocate memory for key */
2763 pgl_keyptr key = palloc(sizeof(pgl_pointkey));
2764 /* transform point to key */
2765 pgl_point_to_key(point, key);
2766 /* create new GISTENTRY structure as return value */
2767 retval = palloc(sizeof(GISTENTRY));
2768 gistentryinit(
2769 *retval, PointerGetDatum(key),
2770 entry->rel, entry->page, entry->offset, FALSE
2771 );
2772 } else {
2773 /* inner nodes have already been transformed */
2774 retval = entry;
2775 }
2776 /* return pointer to old or new GISTENTRY structure */
2777 PG_RETURN_POINTER(retval);
2778 }
2780 /* GiST "compress" support function for indicis on circles */
2781 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
2782 Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
2783 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2784 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2785 /* only transform new leaves */
2786 if (entry->leafkey) {
2787 /* get circle to be transformed */
2788 pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
2789 /* allocate memory for key */
2790 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2791 /* transform circle to key */
2792 pgl_circle_to_key(circle, key);
2793 /* create new GISTENTRY structure as return value */
2794 retval = palloc(sizeof(GISTENTRY));
2795 gistentryinit(
2796 *retval, PointerGetDatum(key),
2797 entry->rel, entry->page, entry->offset, FALSE
2798 );
2799 } else {
2800 /* inner nodes have already been transformed */
2801 retval = entry;
2802 }
2803 /* return pointer to old or new GISTENTRY structure */
2804 PG_RETURN_POINTER(retval);
2805 }
2807 /* GiST "compress" support function for indices on clusters */
2808 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
2809 Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
2810 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2811 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2812 /* only transform new leaves */
2813 if (entry->leafkey) {
2814 /* get cluster to be transformed (detoasting necessary!) */
2815 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
2816 /* allocate memory for key */
2817 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2818 /* transform cluster to key */
2819 pgl_circle_to_key(&(cluster->bounding), key);
2820 /* create new GISTENTRY structure as return value */
2821 retval = palloc(sizeof(GISTENTRY));
2822 gistentryinit(
2823 *retval, PointerGetDatum(key),
2824 entry->rel, entry->page, entry->offset, FALSE
2825 );
2826 /* free detoasted datum */
2827 if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
2828 } else {
2829 /* inner nodes have already been transformed */
2830 retval = entry;
2831 }
2832 /* return pointer to old or new GISTENTRY structure */
2833 PG_RETURN_POINTER(retval);
2834 }
2836 /* GiST "decompress" support function for indices */
2837 PG_FUNCTION_INFO_V1(pgl_gist_decompress);
2838 Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
2839 /* return passed pointer without transformation */
2840 PG_RETURN_POINTER(PG_GETARG_POINTER(0));
2841 }
2843 /* GiST "penalty" support function */
2844 PG_FUNCTION_INFO_V1(pgl_gist_penalty);
2845 Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
2846 GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
2847 GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
2848 float *penalty = (float *)PG_GETARG_POINTER(2);
2849 /* get original key and key to insert */
2850 pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
2851 pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
2852 /* copy original key */
2853 union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
2854 if (PGL_KEY_IS_AREAKEY(orig)) {
2855 memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
2856 } else {
2857 memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
2858 }
2859 /* calculate union of both keys */
2860 pgl_unite_keys((pgl_keyptr)&union_key, new);
2861 /* penalty equal to reduction of key length (logarithm of added area) */
2862 /* (return value by setting referenced value and returning pointer) */
2863 *penalty = (
2864 PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
2865 );
2866 PG_RETURN_POINTER(penalty);
2867 }
2869 /* GiST "picksplit" support function */
2870 PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
2871 Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
2872 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2873 GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
2874 OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */
2875 union {
2876 pgl_pointkey pointkey;
2877 pgl_areakey areakey;
2878 } union_all; /* union of all keys (to be calculated from scratch)
2879 (later cut in half) */
2880 int is_areakey = PGL_KEY_IS_AREAKEY(
2881 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
2882 );
2883 int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
2884 pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */
2885 pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */
2886 pgl_keyptr key; /* current key to be processed */
2887 /* allocate memory for array of left and right keys, set counts to zero */
2888 v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2889 v->spl_nleft = 0;
2890 v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2891 v->spl_nright = 0;
2892 /* calculate union of all keys from scratch */
2893 memcpy(
2894 (pgl_keyptr)&union_all,
2895 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
2896 keysize
2897 );
2898 for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
2899 pgl_unite_keys(
2900 (pgl_keyptr)&union_all,
2901 (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
2902 );
2903 }
2904 /* check if trivial split is necessary due to exhausted key length */
2905 /* (Note: keys for empty objects must have node depth set to maximum) */
2906 if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
2907 is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
2908 )) {
2909 /* half of all keys go left */
2910 for (
2911 i=FirstOffsetNumber;
2912 i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2913 i=OffsetNumberNext(i)
2914 ) {
2915 /* pointer to current key */
2916 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2917 /* update unionL */
2918 /* check if key is first key that goes left */
2919 if (!v->spl_nleft) {
2920 /* first key that goes left is just copied to unionL */
2921 memcpy(unionL, key, keysize);
2922 } else {
2923 /* unite current value and next key */
2924 pgl_unite_keys(unionL, key);
2925 }
2926 /* append offset number to list of keys that go left */
2927 v->spl_left[v->spl_nleft++] = i;
2928 }
2929 /* other half goes right */
2930 for (
2931 i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2932 i<entryvec->n;
2933 i=OffsetNumberNext(i)
2934 ) {
2935 /* pointer to current key */
2936 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2937 /* update unionR */
2938 /* check if key is first key that goes right */
2939 if (!v->spl_nright) {
2940 /* first key that goes right is just copied to unionR */
2941 memcpy(unionR, key, keysize);
2942 } else {
2943 /* unite current value and next key */
2944 pgl_unite_keys(unionR, key);
2945 }
2946 /* append offset number to list of keys that go right */
2947 v->spl_right[v->spl_nright++] = i;
2948 }
2949 }
2950 /* otherwise, a non-trivial split is possible */
2951 else {
2952 /* cut covered area in half */
2953 /* (union_all then refers to area of keys that go left) */
2954 /* check if union of all keys covers empty and non-empty objects */
2955 if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
2956 /* if yes, split into empty and non-empty objects */
2957 pgl_key_set_empty((pgl_keyptr)&union_all);
2958 } else {
2959 /* otherwise split by next bit */
2960 ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
2961 /* NOTE: type bit conserved */
2962 }
2963 /* determine for each key if it goes left or right */
2964 for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
2965 /* pointer to current key */
2966 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2967 /* keys within one half of the area go left */
2968 if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
2969 /* update unionL */
2970 /* check if key is first key that goes left */
2971 if (!v->spl_nleft) {
2972 /* first key that goes left is just copied to unionL */
2973 memcpy(unionL, key, keysize);
2974 } else {
2975 /* unite current value of unionL and processed key */
2976 pgl_unite_keys(unionL, key);
2977 }
2978 /* append offset number to list of keys that go left */
2979 v->spl_left[v->spl_nleft++] = i;
2980 }
2981 /* the other keys go right */
2982 else {
2983 /* update unionR */
2984 /* check if key is first key that goes right */
2985 if (!v->spl_nright) {
2986 /* first key that goes right is just copied to unionR */
2987 memcpy(unionR, key, keysize);
2988 } else {
2989 /* unite current value of unionR and processed key */
2990 pgl_unite_keys(unionR, key);
2991 }
2992 /* append offset number to list of keys that go right */
2993 v->spl_right[v->spl_nright++] = i;
2994 }
2995 }
2996 }
2997 /* store unions in return value */
2998 v->spl_ldatum = PointerGetDatum(unionL);
2999 v->spl_rdatum = PointerGetDatum(unionR);
3000 /* return all results */
3001 PG_RETURN_POINTER(v);
3002 }
3004 /* GiST "same"/"equal" support function */
3005 PG_FUNCTION_INFO_V1(pgl_gist_same);
3006 Datum pgl_gist_same(PG_FUNCTION_ARGS) {
3007 pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
3008 pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
3009 bool *boolptr = (bool *)PG_GETARG_POINTER(2);
3010 /* two keys are equal if they are binary equal */
3011 /* (return result by setting referenced boolean and returning pointer) */
3012 *boolptr = !memcmp(
3013 key1,
3014 key2,
3015 PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
3016 );
3017 PG_RETURN_POINTER(boolptr);
3018 }
3020 /* GiST "distance" support function */
3021 PG_FUNCTION_INFO_V1(pgl_gist_distance);
3022 Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
3023 GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
3024 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
3025 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
3026 bool *recheck = (bool *)PG_GETARG_POINTER(4);
3027 double distance; /* return value */
3028 /* demand recheck because distance is just an estimation */
3029 /* (real distance may be bigger) */
3030 *recheck = true;
3031 /* strategy number aliases for different operators using the same strategy */
3032 strategy %= 100;
3033 /* strategy number 31: distance to point */
3034 if (strategy == 31) {
3035 /* query datum is a point */
3036 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
3037 /* use pgl_estimate_pointkey_distance() function to compute result */
3038 distance = pgl_estimate_key_distance(key, query);
3039 /* avoid infinity (reserved!) */
3040 if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3041 /* return result */
3042 PG_RETURN_FLOAT8(distance);
3043 }
3044 /* strategy number 33: distance to circle */
3045 if (strategy == 33) {
3046 /* query datum is a circle */
3047 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
3048 /* estimate distance to circle center and substract circle radius */
3049 distance = (
3050 pgl_estimate_key_distance(key, &(query->center)) - query->radius
3051 );
3052 /* convert non-positive values to zero and avoid infinity (reserved!) */
3053 if (distance <= 0) distance = 0;
3054 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3055 /* return result */
3056 PG_RETURN_FLOAT8(distance);
3057 }
3058 /* strategy number 34: distance to cluster */
3059 if (strategy == 34) {
3060 /* query datum is a cluster */
3061 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3062 /* estimate distance to bounding center and substract bounding radius */
3063 distance = (
3064 pgl_estimate_key_distance(key, &(query->bounding.center)) -
3065 query->bounding.radius
3066 );
3067 /* convert non-positive values to zero and avoid infinity (reserved!) */
3068 if (distance <= 0) distance = 0;
3069 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
3070 /* free detoasted cluster (if copy) */
3071 PG_FREE_IF_COPY(query, 1);
3072 /* return result */
3073 PG_RETURN_FLOAT8(distance);
3074 }
3075 /* throw error for any unknown strategy number */
3076 elog(ERROR, "unrecognized strategy number: %d", strategy);
3077 }