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