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