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