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