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