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