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