pgLatLon
view latlon-v0001.c @ 2:4f07a22f4d45
Minor change in introduction of README file
| author | jbe |
|---|---|
| date | Sun Aug 21 20:13:05 2016 +0200 (2016-08-21) |
| parents | 3b70e93cc07d |
| children |
line source
2 /*-------------*
3 * C prelude *
4 *-------------*/
6 #include "postgres.h"
7 #include "fmgr.h"
8 #include "libpq/pqformat.h"
9 #include "access/gist.h"
10 #include "access/stratnum.h"
11 #include "utils/array.h"
12 #include <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_DIAMETER ((4.0*PGL_SPHEROID_A + 2.0*PGL_SPHEROID_B) / 3.0)
41 #define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */
42 #define PGL_FADELIMIT (PGL_DIAMETER * M_PI / 6.0) /* 1/6 circumference */
43 #define PGL_MAXDIST (PGL_DIAMETER * M_PI / 2.0) /* maximum distance */
45 /* calculate distance between two points on earth (given in degrees) */
46 static inline double pgl_distance(
47 double lat1, double lon1, double lat2, double lon2
48 ) {
49 float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
50 float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
51 /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
52 /* lon1 = 0 (not used anymore) */
53 lon2 = fabs(lon2-lon1);
54 /* convert to radians (first divide, then multiply) */
55 lat1 = (lat1 / 180.0) * M_PI;
56 lat2 = (lat2 / 180.0) * M_PI;
57 lon2 = (lon2 / 180.0) * M_PI;
58 /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
59 operations (lon2 >= lon1 is already ensured in a previous step) */
60 if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
61 /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
62 of 1.0 */
63 lat1cos = cos(lat1); lat1sin = sin(lat1);
64 lat2cos = cos(lat2); lat2sin = sin(lat2);
65 lon2cos = cos(lon2); lon2sin = sin(lon2);
66 nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
67 nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
68 x1 = nphi1 * lat1cos;
69 z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
70 x2 = nphi2 * lat2cos * lon2cos;
71 y2 = nphi2 * lat2cos * lon2sin;
72 z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
73 /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
74 g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
75 /* convert tunnel distance through scaled ellipsoid to approximated surface
76 distance on original ellipsoid */
77 if (g > 1.0) g = 1.0;
78 s = PGL_DIAMETER * asin(g);
79 /* return result only if small enough to be precise (less than 1/3 of
80 maximum possible distance) */
81 if (s <= PGL_FADELIMIT) return s;
82 /* determine antipodal point of second point (i.e. mirror second point) */
83 lat2 = -lat2; lon2 = lon2 - M_PI;
84 lat2cos = cos(lat2); lat2sin = sin(lat2);
85 lon2cos = cos(lon2); lon2sin = sin(lon2);
86 /* calculate 3d coordinates of antipodal point on scaled ellipsoid */
87 nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
88 x2 = nphi2 * lat2cos * lon2cos;
89 y2 = nphi2 * lat2cos * lon2sin;
90 z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
91 /* calculate tunnel distance to antipodal point through scaled ellipsoid */
92 g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
93 /* convert tunnel distance to antipodal point through scaled ellipsoid to
94 approximated surface distance to antipodal point on original ellipsoid */
95 if (g > 1.0) g = 1.0;
96 t = PGL_DIAMETER * asin(g);
97 /* surface distance between original points can now be approximated by
98 substracting antipodal distance from maximum possible distance;
99 return result only if small enough (less than 1/3 of maximum possible
100 distance) */
101 if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
102 /* otherwise crossfade direct and antipodal result to ensure monotonicity */
103 return (
104 (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
105 (s + t - 2*PGL_FADELIMIT)
106 );
107 }
109 /* finite distance that can not be reached on earth */
110 #define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
113 /*--------------------------------*
114 * simple geographic data types *
115 *--------------------------------*/
117 /* point on earth given by latitude and longitude in degrees */
118 /* (type "epoint" in SQL) */
119 typedef struct {
120 double lat; /* between -90 and 90 (both inclusive) */
121 double lon; /* between -180 and 180 (both inclusive) */
122 } pgl_point;
124 /* box delimited by two parallels and two meridians (all in degrees) */
125 /* (type "ebox" in SQL) */
126 typedef struct {
127 double lat_min; /* between -90 and 90 (both inclusive) */
128 double lat_max; /* between -90 and 90 (both inclusive) */
129 double lon_min; /* between -180 and 180 (both inclusive) */
130 double lon_max; /* between -180 and 180 (both inclusive) */
131 /* if lat_min > lat_max, then box is empty */
132 /* if lon_min > lon_max, then 180th meridian is crossed */
133 } pgl_box;
135 /* circle on earth surface (for radial searches with fixed radius) */
136 /* (type "ecircle" in SQL) */
137 typedef struct {
138 pgl_point center;
139 double radius; /* positive (including +0 but excluding -0), or -INFINITY */
140 /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
141 zero radius (0) denotes a single point,
142 a finite radius (0 < radius < INFINITY) denotes a filled circle, and
143 a radius of INFINITY is valid and means complete coverage of earth. */
144 } pgl_circle;
147 /*----------------------------------*
148 * geographic "cluster" data type *
149 *----------------------------------*/
151 /* A cluster is a collection of points, paths, outlines, and polygons. If two
152 polygons in a cluster overlap, the area covered by both polygons does not
153 belong to the cluster. This way, a cluster can be used to describe complex
154 shapes like polygons with holes. Outlines are non-filled polygons. Paths are
155 open by default (i.e. the last point in the list is not connected with the
156 first point in the list). Note that each outline or polygon in a cluster
157 must cover a longitude range of less than 180 degrees to avoid ambiguities.
158 Areas which are larger may be split into multiple polygons. */
160 /* maximum number of points in a cluster */
161 /* (limited to avoid integer overflows, e.g. when allocating memory) */
162 #define PGL_CLUSTER_MAXPOINTS 16777216
164 /* types of cluster entries */
165 #define PGL_ENTRY_POINT 1 /* a point */
166 #define PGL_ENTRY_PATH 2 /* a path from first point to last point */
167 #define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */
168 #define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */
170 /* Entries of a cluster are described by two different structs: pgl_newentry
171 and pgl_entry. The first is used only during construction of a cluster, the
172 second is used in all other cases (e.g. when reading clusters from the
173 database, performing operations, etc). */
175 /* entry for new geographic cluster during construction of that cluster */
176 typedef struct {
177 int32_t entrytype;
178 int32_t npoints;
179 pgl_point *points; /* pointer to an array of points (pgl_point) */
180 } pgl_newentry;
182 /* entry of geographic cluster */
183 typedef struct {
184 int32_t entrytype; /* type of entry: point, path, outline, polygon */
185 int32_t npoints; /* number of stored points (set to 1 for point entry) */
186 int32_t offset; /* offset of pgl_point array from cluster base address */
187 /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
188 } pgl_entry;
190 /* geographic cluster which is a collection of points, (open) paths, polygons,
191 and outlines (non-filled polygons) */
192 typedef struct {
193 char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */
194 int32_t nentries; /* number of stored points */
195 pgl_circle bounding; /* bounding circle */
196 /* Note: bounding circle ensures alignment of pgl_cluster for points */
197 pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */
198 } pgl_cluster;
200 /* macro to determine memory alignment of points */
201 /* (needed to store pgl_point array after entries in pgl_cluster) */
202 typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
203 #define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
205 /* macro to extract a pointer to the array of points of a cluster entry */
206 #define PGL_ENTRY_POINTS(cluster, idx) \
207 ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
209 /* convert pgl_newentry array to pgl_cluster */
210 static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
211 int i; /* index of current entry */
212 int npoints = 0; /* number of points in whole cluster */
213 int entry_npoints; /* number of points in current entry */
214 int points_offset = PGL_POINT_ALIGNMENT * (
215 ( offsetof(pgl_cluster, entries) +
216 nentries * sizeof(pgl_entry) +
217 PGL_POINT_ALIGNMENT - 1
218 ) / PGL_POINT_ALIGNMENT
219 ); /* offset of pgl_point array from base address (considering alignment) */
220 pgl_cluster *cluster; /* new cluster to be returned */
221 /* determine total number of points */
222 for (i=0; i<nentries; i++) npoints += entries[i].npoints;
223 /* allocate memory for cluster (including entries and points) */
224 cluster = palloc(points_offset + npoints * sizeof(pgl_point));
225 /* re-count total number of points to determine offset for each entry */
226 npoints = 0;
227 /* copy entries and points */
228 for (i=0; i<nentries; i++) {
229 /* determine number of points in entry */
230 entry_npoints = entries[i].npoints;
231 /* copy entry */
232 cluster->entries[i].entrytype = entries[i].entrytype;
233 cluster->entries[i].npoints = entry_npoints;
234 /* calculate offset (in bytes) of pgl_point array */
235 cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
236 /* copy points */
237 memcpy(
238 PGL_ENTRY_POINTS(cluster, i),
239 entries[i].points,
240 entry_npoints * sizeof(pgl_point)
241 );
242 /* update total number of points processed */
243 npoints += entry_npoints;
244 }
245 /* set number of entries in cluster */
246 cluster->nentries = nentries;
247 /* set PostgreSQL header for variable sized data */
248 SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
249 /* return newly created cluster */
250 return cluster;
251 }
254 /*----------------------------------------*
255 * C functions on geographic data types *
256 *----------------------------------------*/
258 /* round latitude or longitude to 12 digits after decimal point */
259 static inline double pgl_round(double val) {
260 return round(val * 1e12) / 1e12;
261 }
263 /* compare two points */
264 /* (equality when same point on earth is described, otherwise an arbitrary
265 linear order) */
266 static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
267 double lon1, lon2; /* modified longitudes for special cases */
268 /* use latitude as first ordering criterion */
269 if (point1->lat < point2->lat) return -1;
270 if (point1->lat > point2->lat) return 1;
271 /* determine modified longitudes (considering special case of poles and
272 180th meridian which can be described as W180 or E180) */
273 if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
274 else if (point1->lon == 180) lon1 = -180;
275 else lon1 = point1->lon;
276 if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
277 else if (point2->lon == 180) lon2 = -180;
278 else lon2 = point2->lon;
279 /* use (modified) longitude as secondary ordering criterion */
280 if (lon1 < lon2) return -1;
281 if (lon1 > lon2) return 1;
282 /* no difference found, points are equal */
283 return 0;
284 }
286 /* compare two boxes */
287 /* (equality when same box on earth is described, otherwise an arbitrary linear
288 order) */
289 static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
290 /* two empty boxes are equal, and an empty box is always considered "less
291 than" a non-empty box */
292 if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
293 if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
294 if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
295 /* use southern border as first ordering criterion */
296 if (box1->lat_min < box2->lat_min) return -1;
297 if (box1->lat_min > box2->lat_min) return 1;
298 /* use northern border as second ordering criterion */
299 if (box1->lat_max < box2->lat_max) return -1;
300 if (box1->lat_max > box2->lat_max) return 1;
301 /* use western border as third ordering criterion */
302 if (box1->lon_min < box2->lon_min) return -1;
303 if (box1->lon_min > box2->lon_min) return 1;
304 /* use eastern border as fourth ordering criterion */
305 if (box1->lon_max < box2->lon_max) return -1;
306 if (box1->lon_max > box2->lon_max) return 1;
307 /* no difference found, boxes are equal */
308 return 0;
309 }
311 /* compare two circles */
312 /* (equality when same circle on earth is described, otherwise an arbitrary
313 linear order) */
314 static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
315 /* two circles with same infinite radius (positive or negative infinity) are
316 considered equal independently of center point */
317 if (
318 !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
319 circle1->radius == circle2->radius
320 ) return 0;
321 /* use radius as first ordering criterion */
322 if (circle1->radius < circle2->radius) return -1;
323 if (circle1->radius > circle2->radius) return 1;
324 /* use center point as secondary ordering criterion */
325 return pgl_point_cmp(&(circle1->center), &(circle2->center));
326 }
328 /* set box to empty box*/
329 static void pgl_box_set_empty(pgl_box *box) {
330 box->lat_min = INFINITY;
331 box->lat_max = -INFINITY;
332 box->lon_min = 0;
333 box->lon_max = 0;
334 }
336 /* check if point is inside a box */
337 static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
338 return (
339 point->lat >= box->lat_min && point->lat <= box->lat_max && (
340 (box->lon_min > box->lon_max) ? (
341 /* box crosses 180th meridian */
342 point->lon >= box->lon_min || point->lon <= box->lon_max
343 ) : (
344 /* box does not cross the 180th meridian */
345 point->lon >= box->lon_min && point->lon <= box->lon_max
346 )
347 )
348 );
349 }
351 /* check if two boxes overlap */
352 static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
353 return (
354 box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */
355 ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
356 ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
357 (
358 /* check if one and only one box crosses the 180th meridian */
359 ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
360 ((box2->lon_min > box2->lon_max) ? 1 : 0)
361 ) ? (
362 /* exactly one box crosses the 180th meridian */
363 box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
364 box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
365 ) : (
366 /* no box or both boxes cross the 180th meridian */
367 (
368 (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
369 (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
370 ) ||
371 /* handle W180 == E180 */
372 ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
373 ( box2->lon_min == -180 && box1->lon_max == 180 )
374 )
375 )
376 );
377 }
379 /* check unambiguousness of east/west orientation of cluster entries and set
380 bounding circle of cluster */
381 static bool pgl_finalize_cluster(pgl_cluster *cluster) {
382 int i, j; /* i: index of entry, j: index of point in entry */
383 int npoints; /* number of points in entry */
384 int total_npoints = 0; /* total number of points in cluster */
385 pgl_point *points; /* points in entry */
386 int lon_dir; /* first point of entry west (-1) or east (+1) */
387 double lon_break = 0; /* antipodal longitude of first point in entry */
388 double lon_min, lon_max; /* covered longitude range of entry */
389 double value; /* temporary variable */
390 /* reset bounding circle center to empty circle at 0/0 coordinates */
391 cluster->bounding.center.lat = 0;
392 cluster->bounding.center.lon = 0;
393 cluster->bounding.radius = -INFINITY;
394 /* if cluster is not empty */
395 if (cluster->nentries != 0) {
396 /* iterate over all cluster entries and ensure they each cover a longitude
397 range less than 180 degrees */
398 for (i=0; i<cluster->nentries; i++) {
399 /* get properties of entry */
400 npoints = cluster->entries[i].npoints;
401 points = PGL_ENTRY_POINTS(cluster, i);
402 /* get longitude of first point of entry */
403 value = points[0].lon;
404 /* initialize lon_min and lon_max with longitude of first point */
405 lon_min = value;
406 lon_max = value;
407 /* determine east/west orientation of first point and calculate antipodal
408 longitude (Note: rounding required here) */
409 if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
410 else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); }
411 else lon_dir = 0;
412 /* iterate over all other points in entry */
413 for (j=1; j<npoints; j++) {
414 /* consider longitude wrap-around */
415 value = points[j].lon;
416 if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
417 else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
418 /* update lon_min and lon_max */
419 if (value < lon_min) lon_min = value;
420 else if (value > lon_max) lon_max = value;
421 /* return false if 180 degrees or more are covered */
422 if (lon_max - lon_min >= 180) return false;
423 }
424 }
425 /* iterate over all points of all entries and calculate arbitrary center
426 point for bounding circle (best if center point minimizes the radius,
427 but some error is allowed here) */
428 for (i=0; i<cluster->nentries; i++) {
429 /* get properties of entry */
430 npoints = cluster->entries[i].npoints;
431 points = PGL_ENTRY_POINTS(cluster, i);
432 /* check if first entry */
433 if (i==0) {
434 /* get longitude of first point of first entry in whole cluster */
435 value = points[0].lon;
436 /* initialize lon_min and lon_max with longitude of first point of
437 first entry in whole cluster (used to determine if whole cluster
438 covers a longitude range of 180 degrees or more) */
439 lon_min = value;
440 lon_max = value;
441 /* determine east/west orientation of first point and calculate
442 antipodal longitude (Note: rounding not necessary here) */
443 if (value < 0) { lon_dir = -1; lon_break = value + 180; }
444 else if (value > 0) { lon_dir = 1; lon_break = value - 180; }
445 else lon_dir = 0;
446 }
447 /* iterate over all points in entry */
448 for (j=0; j<npoints; j++) {
449 /* longitude wrap-around (Note: rounding not necessary here) */
450 value = points[j].lon;
451 if (lon_dir < 0 && value > lon_break) value -= 360;
452 else if (lon_dir > 0 && value < lon_break) value += 360;
453 if (value < lon_min) lon_min = value;
454 else if (value > lon_max) lon_max = value;
455 /* set bounding circle to cover whole earth if more than 180 degrees
456 are covered */
457 if (lon_max - lon_min >= 180) {
458 cluster->bounding.center.lat = 0;
459 cluster->bounding.center.lon = 0;
460 cluster->bounding.radius = INFINITY;
461 return true;
462 }
463 /* add point to bounding circle center (for average calculation) */
464 cluster->bounding.center.lat += points[j].lat;
465 cluster->bounding.center.lon += value;
466 }
467 /* count total number of points */
468 total_npoints += npoints;
469 }
470 /* determine average latitude and longitude of cluster */
471 cluster->bounding.center.lat /= total_npoints;
472 cluster->bounding.center.lon /= total_npoints;
473 /* normalize longitude of center of cluster bounding circle */
474 if (cluster->bounding.center.lon < -180) {
475 cluster->bounding.center.lon += 360;
476 }
477 else if (cluster->bounding.center.lon > 180) {
478 cluster->bounding.center.lon -= 360;
479 }
480 /* round bounding circle center (useful if it is used by other functions) */
481 cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
482 cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
483 /* calculate radius of bounding circle */
484 for (i=0; i<cluster->nentries; i++) {
485 npoints = cluster->entries[i].npoints;
486 points = PGL_ENTRY_POINTS(cluster, i);
487 for (j=0; j<npoints; j++) {
488 value = pgl_distance(
489 cluster->bounding.center.lat, cluster->bounding.center.lon,
490 points[j].lat, points[j].lon
491 );
492 if (value > cluster->bounding.radius) cluster->bounding.radius = value;
493 }
494 }
495 }
496 /* return true (east/west orientation is unambiguous) */
497 return true;
498 }
500 /* check if point is inside cluster */
501 static bool pgl_point_in_cluster(pgl_point *point, pgl_cluster *cluster) {
502 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
503 int entrytype; /* type of entry */
504 int npoints; /* number of points in entry */
505 pgl_point *points; /* array of points in entry */
506 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
507 double lon_break = 0; /* antipodal longitude of first vertex */
508 double lat0 = point->lat; /* latitude of point */
509 double lon0; /* (adjusted) longitude of point */
510 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
511 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
512 double lon; /* longitude of intersection */
513 int counter = 0; /* counter for intersections east of point */
514 /* points outside bounding circle are always assumed to be non-overlapping */
515 /* (necessary for consistent table and index scans) */
516 if (
517 pgl_distance(
518 point->lat, point->lon,
519 cluster->bounding.center.lat, cluster->bounding.center.lon
520 ) > cluster->bounding.radius
521 ) return false;
522 /* iterate over all entries */
523 for (i=0; i<cluster->nentries; i++) {
524 /* get properties of entry */
525 entrytype = cluster->entries[i].entrytype;
526 npoints = cluster->entries[i].npoints;
527 points = PGL_ENTRY_POINTS(cluster, i);
528 /* determine east/west orientation of first point of entry and calculate
529 antipodal longitude */
530 lon_break = points[0].lon;
531 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
532 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
533 else lon_dir = 0;
534 /* get longitude of point */
535 lon0 = point->lon;
536 /* consider longitude wrap-around for point */
537 if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
538 else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
539 /* iterate over all edges and vertices */
540 for (j=0; j<npoints; j++) {
541 /* return true if point is on vertex of polygon */
542 if (pgl_point_cmp(point, &(points[j])) == 0) return true;
543 /* calculate index of next vertex */
544 k = (j+1) % npoints;
545 /* skip last edge unless entry is (closed) outline or polygon */
546 if (
547 k == 0 &&
548 entrytype != PGL_ENTRY_OUTLINE &&
549 entrytype != PGL_ENTRY_POLYGON
550 ) continue;
551 /* get latitude and longitude values of edge */
552 lat1 = points[j].lat;
553 lat2 = points[k].lat;
554 lon1 = points[j].lon;
555 lon2 = points[k].lon;
556 /* consider longitude wrap-around for edge */
557 if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
558 else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
559 if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
560 else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
561 /* return true if point is on horizontal (west to east) edge of polygon */
562 if (
563 lat0 == lat1 && lat0 == lat2 &&
564 ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
565 ) return true;
566 /* check if edge crosses east/west line of point */
567 if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
568 /* calculate longitude of intersection */
569 lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
570 /* return true if intersection goes (approximately) through point */
571 if (pgl_round(lon) == lon0) return true;
572 /* count intersection if east of point and entry is polygon*/
573 if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
574 }
575 }
576 }
577 /* return true if number of intersections is odd */
578 return counter & 1;
579 }
581 /* calculate (approximate) distance between point and cluster */
582 static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
583 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
584 int entrytype; /* type of entry */
585 int npoints; /* number of points in entry */
586 pgl_point *points; /* array of points in entry */
587 int lon_dir = 0; /* first vertex west (-1) or east (+1) */
588 double lon_break = 0; /* antipodal longitude of first vertex */
589 double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */
590 double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */
591 double lat0 = point->lat; /* latitude of point */
592 double lon0; /* (adjusted) longitude of point */
593 double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
594 double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
595 double s; /* scalar for vector calculations */
596 double dist; /* distance calculated in one step */
597 double min_dist = INFINITY; /* minimum distance */
598 /* distance is zero if point is contained in cluster */
599 if (pgl_point_in_cluster(point, cluster)) return 0;
600 /* iterate over all entries */
601 for (i=0; i<cluster->nentries; i++) {
602 /* get properties of entry */
603 entrytype = cluster->entries[i].entrytype;
604 npoints = cluster->entries[i].npoints;
605 points = PGL_ENTRY_POINTS(cluster, i);
606 /* determine east/west orientation of first point of entry and calculate
607 antipodal longitude */
608 lon_break = points[0].lon;
609 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
610 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
611 else lon_dir = 0;
612 /* determine covered longitude range */
613 for (j=0; j<npoints; j++) {
614 /* get longitude of vertex */
615 lon1 = points[j].lon;
616 /* adjust longitude to fix potential wrap-around */
617 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
618 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
619 /* update minimum and maximum longitude of polygon */
620 if (j == 0 || lon1 < lon_min) lon_min = lon1;
621 if (j == 0 || lon1 > lon_max) lon_max = lon1;
622 }
623 /* adjust longitude wrap-around according to full longitude range */
624 lon_break = (lon_max + lon_min) / 2;
625 if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
626 else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
627 /* get longitude of point */
628 lon0 = point->lon;
629 /* consider longitude wrap-around for point */
630 if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
631 else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
632 /* iterate over all edges and vertices */
633 for (j=0; j<npoints; j++) {
634 /* get latitude and longitude values of current point */
635 lat1 = points[j].lat;
636 lon1 = points[j].lon;
637 /* consider longitude wrap-around for current point */
638 if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
639 else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
640 /* calculate distance to vertex */
641 dist = pgl_distance(lat0, lon0, lat1, lon1);
642 /* store calculated distance if smallest */
643 if (dist < min_dist) min_dist = dist;
644 /* calculate index of next vertex */
645 k = (j+1) % npoints;
646 /* skip last edge unless entry is (closed) outline or polygon */
647 if (
648 k == 0 &&
649 entrytype != PGL_ENTRY_OUTLINE &&
650 entrytype != PGL_ENTRY_POLYGON
651 ) continue;
652 /* get latitude and longitude values of next point */
653 lat2 = points[k].lat;
654 lon2 = points[k].lon;
655 /* consider longitude wrap-around for next point */
656 if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
657 else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
658 /* go to next vertex and edge if edge is degenerated */
659 if (lat1 == lat2 && lon1 == lon2) continue;
660 /* otherwise test if point can be projected onto edge of polygon */
661 s = (
662 ((lat0-lat1) * (lat2-lat1) + (lon0-lon1) * (lon2-lon1)) /
663 ((lat2-lat1) * (lat2-lat1) + (lon2-lon1) * (lon2-lon1))
664 );
665 /* go to next vertex and edge if point cannot be projected */
666 if (!(s > 0 && s < 1)) continue;
667 /* calculate distance from original point to projected point */
668 dist = pgl_distance(
669 lat0, lon0,
670 lat1 + s * (lat2-lat1),
671 lon1 + s * (lon2-lon1)
672 );
673 /* store calculated distance if smallest */
674 if (dist < min_dist) min_dist = dist;
675 }
676 }
677 /* return minimum distance */
678 return min_dist;
679 }
681 /* estimator function for distance between box and point */
682 /* allowed to return smaller values than actually correct */
683 static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
684 double dlon; /* longitude range of box (delta longitude) */
685 double h; /* half of distance along meridian */
686 double d; /* distance between both southern or both northern points */
687 double cur_dist; /* calculated distance */
688 double min_dist; /* minimum distance calculated */
689 /* return infinity if bounding box is empty */
690 if (box->lat_min > box->lat_max) return INFINITY;
691 /* return zero if point is inside bounding box */
692 if (pgl_point_in_box(point, box)) return 0;
693 /* calculate delta longitude */
694 dlon = box->lon_max - box->lon_min;
695 if (dlon < 0) dlon += 360; /* 180th meridian crossed */
696 /* if delta longitude is greater than 180 degrees, perform safe fall-back */
697 if (dlon > 180) return 0;
698 /* calculate half of distance along meridian */
699 h = pgl_distance(box->lat_min, 0, box->lat_max, 0) / 2;
700 /* calculate full distance between southern points */
701 d = pgl_distance(box->lat_min, 0, box->lat_min, dlon);
702 /* calculate maximum of full distance and half distance */
703 if (h > d) d = h;
704 /* calculate distance from point to first southern vertex and substract
705 maximum error */
706 min_dist = pgl_distance(
707 point->lat, point->lon, box->lat_min, box->lon_min
708 ) - d;
709 /* return zero if estimated distance is smaller than zero */
710 if (min_dist <= 0) return 0;
711 /* repeat procedure with second southern vertex */
712 cur_dist = pgl_distance(
713 point->lat, point->lon, box->lat_min, box->lon_max
714 ) - d;
715 if (cur_dist <= 0) return 0;
716 if (cur_dist < min_dist) min_dist = cur_dist;
717 /* calculate full distance between northern points */
718 d = pgl_distance(box->lat_max, 0, box->lat_max, dlon);
719 /* calculate maximum of full distance and half distance */
720 if (h > d) d = h;
721 /* repeat procedure with northern vertices */
722 cur_dist = pgl_distance(
723 point->lat, point->lon, box->lat_max, box->lon_max
724 ) - d;
725 if (cur_dist <= 0) return 0;
726 if (cur_dist < min_dist) min_dist = cur_dist;
727 cur_dist = pgl_distance(
728 point->lat, point->lon, box->lat_max, box->lon_min
729 ) - d;
730 if (cur_dist <= 0) return 0;
731 if (cur_dist < min_dist) min_dist = cur_dist;
732 /* return smallest value (unless already returned zero) */
733 return min_dist;
734 }
737 /*----------------------------*
738 * fractal geographic index *
739 *----------------------------*/
741 /* number of bytes used for geographic (center) position in keys */
742 #define PGL_KEY_LATLON_BYTELEN 7
744 /* maximum reference value for logarithmic size of geographic objects */
745 #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */
747 /* safety margin to avoid floating point errors in distance estimation */
748 #define PGL_FPE_SAFETY (1.0+1e-14) /* slightly greater than 1.0 */
750 /* pointer to index key (either pgl_pointkey or pgl_areakey) */
751 typedef unsigned char *pgl_keyptr;
753 /* index key for points (objects with zero area) on the spheroid */
754 /* bit 0..55: interspersed bits of latitude and longitude,
755 bit 56..57: always zero,
756 bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
757 typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
759 /* index key for geographic objects on spheroid with area greater than zero */
760 /* bit 0..55: interspersed bits of latitude and longitude of center point,
761 bit 56: always set to 1,
762 bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
763 bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
764 PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
765 = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
766 (with interspersed bits = 0 and node depth = 0) for keys which
767 cover both empty and non-empty objects */
769 typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
771 /* helper macros for reading/writing index keys */
772 #define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN
773 #define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1)
774 #define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8)
775 #define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1)
776 #define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1)
777 #define PGL_AREAKEY_TYPEMASK 0x80
778 #define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
779 #define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
780 ( PGL_KEY_LATLONBIT(key1, n) ^ \
781 PGL_KEY_LATLONBIT(key2, n) )
782 #define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
783 PGL_AREAKEY_TYPEMASK)
784 #define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
785 (PGL_AREAKEY_TYPEMASK-1))
786 #define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET])
787 #define PGL_KEY_OBJSIZE_EMPTY 126
788 #define PGL_KEY_OBJSIZE_UNIVERSAL 127
789 #define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \
790 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
791 PGL_KEY_OBJSIZE_EMPTY )
792 #define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
793 (key)[PGL_KEY_OBJSIZE_OFFSET] == \
794 PGL_KEY_OBJSIZE_UNIVERSAL )
796 /* set area key to match empty objects only */
797 static void pgl_key_set_empty(pgl_keyptr key) {
798 memset(key, 0, sizeof(pgl_areakey));
799 /* Note: setting node depth to maximum is required for picksplit function */
800 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
801 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
802 }
804 /* set area key to match any object (including empty objects) */
805 static void pgl_key_set_universal(pgl_keyptr key) {
806 memset(key, 0, sizeof(pgl_areakey));
807 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
808 key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
809 }
811 /* convert a point on earth into a max-depth key to be used in index */
812 static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
813 double lat = point->lat;
814 double lon = point->lon;
815 int i;
816 /* clear latitude and longitude bits */
817 memset(key, 0, PGL_KEY_LATLON_BYTELEN);
818 /* set node depth to maximum and type bit to zero */
819 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
820 /* iterate over all latitude/longitude bit pairs */
821 for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
822 /* determine latitude bit */
823 if (lat >= 0) {
824 key[i/4] |= 0x80 >> (2*(i%4));
825 lat *= 2; lat -= 90;
826 } else {
827 lat *= 2; lat += 90;
828 }
829 /* determine longitude bit */
830 if (lon >= 0) {
831 key[i/4] |= 0x80 >> (2*(i%4)+1);
832 lon *= 2; lon -= 180;
833 } else {
834 lon *= 2; lon += 180;
835 }
836 }
837 }
839 /* convert a circle on earth into a max-depth key to be used in an index */
840 static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
841 /* handle special case of empty circle */
842 if (circle->radius < 0) {
843 pgl_key_set_empty(key);
844 return;
845 }
846 /* perform same action as for point keys */
847 pgl_point_to_key(&(circle->center), key);
848 /* but overwrite type and node depth to fit area index key */
849 key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
850 /* check if radius is greater than (or equal to) reference size */
851 /* (treat equal values as greater values for numerical safety) */
852 if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
853 /* if yes, set logarithmic size to zero */
854 key[PGL_KEY_OBJSIZE_OFFSET] = 0;
855 } else {
856 /* otherwise, determine logarithmic size iteratively */
857 /* (one step is equivalent to a factor of sqrt(2)) */
858 double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
859 int objsize = 1;
860 while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
861 /* stop when radius is greater than (or equal to) adjusted reference */
862 /* (treat equal values as greater values for numerical safety) */
863 if (circle->radius >= reference) break;
864 reference /= M_SQRT2;
865 objsize++;
866 }
867 /* set logarithmic size to determined value */
868 key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
869 }
870 }
872 /* check if one key is subkey of another key or vice versa */
873 static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
874 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
875 /* determine smallest depth */
876 int depth1 = PGL_KEY_NODEDEPTH(key1);
877 int depth2 = PGL_KEY_NODEDEPTH(key2);
878 int depth = (depth1 < depth2) ? depth1 : depth2;
879 /* check if keys are area keys (assuming that both keys have same type) */
880 if (PGL_KEY_IS_AREAKEY(key1)) {
881 int j = 0; /* bit offset for logarithmic object size bits */
882 int k = 0; /* bit offset for latitude and longitude */
883 /* fetch logarithmic object size information */
884 int objsize1 = PGL_KEY_OBJSIZE(key1);
885 int objsize2 = PGL_KEY_OBJSIZE(key2);
886 /* handle special cases for empty objects (universal and empty keys) */
887 if (
888 objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
889 objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
890 ) return true;
891 if (
892 objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
893 objsize2 == PGL_KEY_OBJSIZE_EMPTY
894 ) return objsize1 == objsize2;
895 /* iterate through key bits */
896 for (i=0; i<depth; i++) {
897 /* every second bit is a bit describing the object size */
898 if (i%2 == 0) {
899 /* check if object size bit is different in both keys (objsize1 and
900 objsize2 describe the minimum index when object size bit is set) */
901 if (
902 (objsize1 <= j && objsize2 > j) ||
903 (objsize2 <= j && objsize1 > j)
904 ) {
905 /* bit differs, therefore keys are in separate branches */
906 return false;
907 }
908 /* increase bit counter for object size bits */
909 j++;
910 }
911 /* all other bits describe latitude and longitude */
912 else {
913 /* check if bit differs in both keys */
914 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
915 /* bit differs, therefore keys are in separate branches */
916 return false;
917 }
918 /* increase bit counter for latitude/longitude bits */
919 k++;
920 }
921 }
922 }
923 /* if not, keys are point keys */
924 else {
925 /* iterate through key bits */
926 for (i=0; i<depth; i++) {
927 /* check if bit differs in both keys */
928 if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
929 /* bit differs, therefore keys are in separate branches */
930 return false;
931 }
932 }
933 }
934 /* return true because keys are in the same branch */
935 return true;
936 }
938 /* combine two keys into new key which covers both original keys */
939 /* (result stored in first argument) */
940 static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
941 int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
942 /* determine smallest depth */
943 int depth1 = PGL_KEY_NODEDEPTH(dst);
944 int depth2 = PGL_KEY_NODEDEPTH(src);
945 int depth = (depth1 < depth2) ? depth1 : depth2;
946 /* check if keys are area keys (assuming that both keys have same type) */
947 if (PGL_KEY_IS_AREAKEY(dst)) {
948 pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */
949 int j = 0; /* bit offset for logarithmic object size bits */
950 int k = 0; /* bit offset for latitude and longitude */
951 /* fetch logarithmic object size information */
952 int objsize1 = PGL_KEY_OBJSIZE(dst);
953 int objsize2 = PGL_KEY_OBJSIZE(src);
954 /* handle special cases for empty objects (universal and empty keys) */
955 if (
956 objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
957 objsize2 > PGL_AREAKEY_MAXOBJSIZE
958 ) {
959 if (
960 objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
961 objsize2 == PGL_KEY_OBJSIZE_EMPTY
962 ) pgl_key_set_empty(dst);
963 else pgl_key_set_universal(dst);
964 return;
965 }
966 /* iterate through key bits */
967 for (i=0; i<depth; i++) {
968 /* every second bit is a bit describing the object size */
969 if (i%2 == 0) {
970 /* increase bit counter for object size bits first */
971 /* (handy when setting objsize variable) */
972 j++;
973 /* check if object size bit is set in neither key */
974 if (objsize1 >= j && objsize2 >= j) {
975 /* set objsize in destination buffer to indicate that size bit is
976 unset in destination buffer at the current bit position */
977 dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
978 }
979 /* break if object size bit is set in one key only */
980 else if (objsize1 >= j || objsize2 >= j) break;
981 }
982 /* all other bits describe latitude and longitude */
983 else {
984 /* break if bit differs in both keys */
985 if (PGL_KEY_LATLONBIT(dst, k)) {
986 if (!PGL_KEY_LATLONBIT(src, k)) break;
987 /* but set bit in destination buffer if bit is set in both keys */
988 dstbuf[k/8] |= 0x80 >> (k%8);
989 } else if (PGL_KEY_LATLONBIT(src, k)) break;
990 /* increase bit counter for latitude/longitude bits */
991 k++;
992 }
993 }
994 /* set common node depth and type bit (type bit = 1) */
995 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
996 /* copy contents of destination buffer to first key */
997 memcpy(dst, dstbuf, sizeof(pgl_areakey));
998 }
999 /* if not, keys are point keys */
1000 else {
1001 pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */
1002 /* iterate through key bits */
1003 for (i=0; i<depth; i++) {
1004 /* break if bit differs in both keys */
1005 if (PGL_KEY_LATLONBIT(dst, i)) {
1006 if (!PGL_KEY_LATLONBIT(src, i)) break;
1007 /* but set bit in destination buffer if bit is set in both keys */
1008 dstbuf[i/8] |= 0x80 >> (i%8);
1009 } else if (PGL_KEY_LATLONBIT(src, i)) break;
1010 }
1011 /* set common node depth (type bit = 0) */
1012 dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
1013 /* copy contents of destination buffer to first key */
1014 memcpy(dst, dstbuf, sizeof(pgl_pointkey));
1015 }
1016 }
1018 /* determine center(!) boundaries and radius estimation of index key */
1019 static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
1020 int i;
1021 /* determine node depth */
1022 int depth = PGL_KEY_NODEDEPTH(key);
1023 /* center point of possible result */
1024 double lat = 0;
1025 double lon = 0;
1026 /* maximum distance of real center point from key center */
1027 double dlat = 90;
1028 double dlon = 180;
1029 /* maximum radius of contained objects */
1030 double radius = 0; /* always return zero for point index keys */
1031 /* check if key is area key */
1032 if (PGL_KEY_IS_AREAKEY(key)) {
1033 /* get logarithmic object size */
1034 int objsize = PGL_KEY_OBJSIZE(key);
1035 /* handle special cases for empty objects (universal and empty keys) */
1036 if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
1037 pgl_box_set_empty(box);
1038 return 0;
1039 } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
1040 box->lat_min = -90;
1041 box->lat_max = 90;
1042 box->lon_min = -180;
1043 box->lon_max = 180;
1044 return 0; /* any value >= 0 would do */
1045 }
1046 /* calculate maximum possible radius of objects covered by the given key */
1047 if (objsize == 0) radius = INFINITY;
1048 else {
1049 radius = PGL_AREAKEY_REFOBJSIZE;
1050 while (--objsize) radius /= M_SQRT2;
1051 }
1052 /* iterate over latitude and longitude bits in key */
1053 /* (every second bit is a latitude or longitude bit) */
1054 for (i=0; i<depth/2; i++) {
1055 /* check if latitude bit */
1056 if (i%2 == 0) {
1057 /* cut latitude dimension in half */
1058 dlat /= 2;
1059 /* increase center latitude if bit is 1, otherwise decrease */
1060 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1061 else lat -= dlat;
1062 }
1063 /* otherwise longitude bit */
1064 else {
1065 /* cut longitude dimension in half */
1066 dlon /= 2;
1067 /* increase center longitude if bit is 1, otherwise decrease */
1068 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1069 else lon -= dlon;
1070 }
1071 }
1072 }
1073 /* if not, keys are point keys */
1074 else {
1075 /* iterate over all bits in key */
1076 for (i=0; i<depth; i++) {
1077 /* check if latitude bit */
1078 if (i%2 == 0) {
1079 /* cut latitude dimension in half */
1080 dlat /= 2;
1081 /* increase center latitude if bit is 1, otherwise decrease */
1082 if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
1083 else lat -= dlat;
1084 }
1085 /* otherwise longitude bit */
1086 else {
1087 /* cut longitude dimension in half */
1088 dlon /= 2;
1089 /* increase center longitude if bit is 1, otherwise decrease */
1090 if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
1091 else lon -= dlon;
1092 }
1093 }
1094 }
1095 /* calculate boundaries from center point and remaining dlat and dlon */
1096 /* (return values through pointer to box) */
1097 box->lat_min = lat - dlat;
1098 box->lat_max = lat + dlat;
1099 box->lon_min = lon - dlon;
1100 box->lon_max = lon + dlon;
1101 /* return radius (as a function return value) */
1102 return radius;
1103 }
1105 /* estimator function for distance between point and index key */
1106 /* allowed to return smaller values than actually correct */
1107 static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
1108 pgl_box box; /* center(!) bounding box of area index key */
1109 /* calculate center(!) bounding box and maximum radius of objects covered
1110 by area index key (radius is zero for point index keys) */
1111 double distance = pgl_key_to_box(key, &box);
1112 /* calculate estimated distance between bounding box of center point of
1113 indexed object and point passed as second argument, then substract maximum
1114 radius of objects covered by index key */
1115 /* (use PGL_FPE_SAFETY factor to cope with minor floating point errors) */
1116 distance = (
1117 pgl_estimate_point_box_distance(point, &box) / PGL_FPE_SAFETY -
1118 distance * PGL_FPE_SAFETY
1119 );
1120 /* truncate negative results to zero */
1121 if (distance <= 0) distance = 0;
1122 /* return result */
1123 return distance;
1124 }
1127 /*---------------------------------*
1128 * helper functions for text I/O *
1129 *---------------------------------*/
1131 #define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */
1133 /* convert floating point number to string (round-trip safe) */
1134 static void pgl_print_float(char *buf, double flt) {
1135 /* check if number is integral */
1136 if (trunc(flt) == flt) {
1137 /* for integral floats use maximum precision */
1138 snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1139 } else {
1140 /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
1141 snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
1142 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
1143 if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
1144 }
1145 }
1147 /* convert latitude floating point number (in degrees) to string */
1148 static void pgl_print_lat(char *buf, double lat) {
1149 if (signbit(lat)) {
1150 /* treat negative latitudes (including -0) as south */
1151 snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
1152 } else {
1153 /* treat positive latitudes (including +0) as north */
1154 snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
1155 }
1156 }
1158 /* convert longitude floating point number (in degrees) to string */
1159 static void pgl_print_lon(char *buf, double lon) {
1160 if (signbit(lon)) {
1161 /* treat negative longitudes (including -0) as west */
1162 snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
1163 } else {
1164 /* treat positive longitudes (including +0) as east */
1165 snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
1166 }
1167 }
1169 /* bit masks used as return value of pgl_scan() function */
1170 #define PGL_SCAN_NONE 0 /* no value has been parsed */
1171 #define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */
1172 #define PGL_SCAN_LON (1<<1) /* longitude has been parsed */
1173 #define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */
1175 /* parse a coordinate (can be latitude or longitude) */
1176 static int pgl_scan(char **str, double *lat, double *lon) {
1177 double val;
1178 int len;
1179 if (
1180 sscanf(*str, " N %lf %n", &val, &len) ||
1181 sscanf(*str, " n %lf %n", &val, &len)
1182 ) {
1183 *str += len; *lat = val; return PGL_SCAN_LAT;
1184 }
1185 if (
1186 sscanf(*str, " S %lf %n", &val, &len) ||
1187 sscanf(*str, " s %lf %n", &val, &len)
1188 ) {
1189 *str += len; *lat = -val; return PGL_SCAN_LAT;
1190 }
1191 if (
1192 sscanf(*str, " E %lf %n", &val, &len) ||
1193 sscanf(*str, " e %lf %n", &val, &len)
1194 ) {
1195 *str += len; *lon = val; return PGL_SCAN_LON;
1196 }
1197 if (
1198 sscanf(*str, " W %lf %n", &val, &len) ||
1199 sscanf(*str, " w %lf %n", &val, &len)
1200 ) {
1201 *str += len; *lon = -val; return PGL_SCAN_LON;
1202 }
1203 return PGL_SCAN_NONE;
1204 }
1207 /*-----------------*
1208 * SQL functions *
1209 *-----------------*/
1211 /* Note: These function names use "epoint", "ebox", etc. notation here instead
1212 of "point", "box", etc. in order to distinguish them from any previously
1213 defined functions. */
1215 /* function needed for dummy types and/or not implemented features */
1216 PG_FUNCTION_INFO_V1(pgl_notimpl);
1217 Datum pgl_notimpl(PG_FUNCTION_ARGS) {
1218 ereport(ERROR, (errmsg("not implemented by pgLatLon")));
1219 }
1221 /* set point to latitude and longitude (including checks) */
1222 static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
1223 /* reject infinite or NaN values */
1224 if (!isfinite(lat) || !isfinite(lon)) {
1225 ereport(ERROR, (
1226 errcode(ERRCODE_DATA_EXCEPTION),
1227 errmsg("epoint requires finite coordinates")
1228 ));
1229 }
1230 /* check latitude bounds */
1231 if (lat < -90) {
1232 ereport(WARNING, (errmsg("latitude exceeds south pole")));
1233 lat = -90;
1234 } else if (lat > 90) {
1235 ereport(WARNING, (errmsg("latitude exceeds north pole")));
1236 lat = 90;
1237 }
1238 /* check longitude bounds */
1239 if (lon < -180) {
1240 ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
1241 lon += 360 - trunc(lon / 360) * 360;
1242 } else if (lon > 180) {
1243 ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
1244 lon -= 360 + trunc(lon / 360) * 360;
1245 }
1246 /* store rounded latitude/longitude values for round-trip safety */
1247 point->lat = pgl_round(lat);
1248 point->lon = pgl_round(lon);
1249 }
1251 /* create point ("epoint" in SQL) from latitude and longitude */
1252 PG_FUNCTION_INFO_V1(pgl_create_epoint);
1253 Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
1254 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
1255 pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
1256 PG_RETURN_POINTER(point);
1257 }
1259 /* parse point ("epoint" in SQL) */
1260 /* format: '[NS]<float> [EW]<float>' */
1261 PG_FUNCTION_INFO_V1(pgl_epoint_in);
1262 Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
1263 char *str = PG_GETARG_CSTRING(0); /* input string */
1264 char *strptr = str; /* current position within string */
1265 int done = 0; /* bit mask storing if latitude or longitude was read */
1266 double lat, lon; /* parsed values as double precision floats */
1267 pgl_point *point; /* return value (to be palloc'ed) */
1268 /* parse two floats (each latitude or longitude) separated by white-space */
1269 done |= pgl_scan(&strptr, &lat, &lon);
1270 if (strptr != str && isspace(strptr[-1])) {
1271 done |= pgl_scan(&strptr, &lat, &lon);
1272 }
1273 /* require end of string, and latitude and longitude parsed successfully */
1274 if (strptr[0] || done != PGL_SCAN_LATLON) {
1275 ereport(ERROR, (
1276 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1277 errmsg("invalid input syntax for type epoint: \"%s\"", str)
1278 ));
1279 }
1280 /* allocate memory for result */
1281 point = (pgl_point *)palloc(sizeof(pgl_point));
1282 /* set latitude and longitude (and perform checks) */
1283 pgl_epoint_set_latlon(point, lat, lon);
1284 /* return result */
1285 PG_RETURN_POINTER(point);
1286 }
1288 /* create box ("ebox" in SQL) that is empty */
1289 PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
1290 Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
1291 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1292 pgl_box_set_empty(box);
1293 PG_RETURN_POINTER(box);
1294 }
1296 /* set box to given boundaries (including checks) */
1297 static void pgl_ebox_set_boundaries(
1298 pgl_box *box,
1299 double lat_min, double lat_max, double lon_min, double lon_max
1300 ) {
1301 /* if minimum latitude is greater than maximum latitude, return empty box */
1302 if (lat_min > lat_max) {
1303 pgl_box_set_empty(box);
1304 return;
1305 }
1306 /* otherwise reject infinite or NaN values */
1307 if (
1308 !isfinite(lat_min) || !isfinite(lat_max) ||
1309 !isfinite(lon_min) || !isfinite(lon_max)
1310 ) {
1311 ereport(ERROR, (
1312 errcode(ERRCODE_DATA_EXCEPTION),
1313 errmsg("ebox requires finite coordinates")
1314 ));
1315 }
1316 /* check latitude bounds */
1317 if (lat_max < -90) {
1318 ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
1319 lat_max = -90;
1320 } else if (lat_max > 90) {
1321 ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
1322 lat_max = 90;
1323 }
1324 if (lat_min < -90) {
1325 ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
1326 lat_min = -90;
1327 } else if (lat_min > 90) {
1328 ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
1329 lat_min = 90;
1330 }
1331 /* check if all longitudes are included */
1332 if (lon_max - lon_min >= 360) {
1333 if (lon_max - lon_min > 360) ereport(WARNING, (
1334 errmsg("longitude coverage greater than 360 degrees")
1335 ));
1336 lon_min = -180;
1337 lon_max = 180;
1338 } else {
1339 /* normalize longitude bounds */
1340 if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
1341 else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360;
1342 if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
1343 else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360;
1344 }
1345 /* store rounded latitude/longitude values for round-trip safety */
1346 box->lat_min = pgl_round(lat_min);
1347 box->lat_max = pgl_round(lat_max);
1348 box->lon_min = pgl_round(lon_min);
1349 box->lon_max = pgl_round(lon_max);
1350 /* ensure that rounding does not change orientation */
1351 if (lon_min > lon_max && box->lon_min == box->lon_max) {
1352 box->lon_min = -180;
1353 box->lon_max = 180;
1354 }
1355 }
1357 /* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
1358 PG_FUNCTION_INFO_V1(pgl_create_ebox);
1359 Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
1360 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1361 pgl_ebox_set_boundaries(
1362 box,
1363 PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
1364 PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
1365 );
1366 PG_RETURN_POINTER(box);
1367 }
1369 /* create box ("ebox" in SQL) from two points ("epoint"s) */
1370 /* (can not be used to cover a longitude range of more than 120 degrees) */
1371 PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
1372 Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
1373 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
1374 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
1375 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1376 double lat_min, lat_max, lon_min, lon_max;
1377 double dlon; /* longitude range (delta longitude) */
1378 /* order latitude and longitude boundaries */
1379 if (point2->lat < point1->lat) {
1380 lat_min = point2->lat;
1381 lat_max = point1->lat;
1382 } else {
1383 lat_min = point1->lat;
1384 lat_max = point2->lat;
1385 }
1386 if (point2->lon < point1->lon) {
1387 lon_min = point2->lon;
1388 lon_max = point1->lon;
1389 } else {
1390 lon_min = point1->lon;
1391 lon_max = point2->lon;
1392 }
1393 /* calculate longitude range (round to avoid floating point errors) */
1394 dlon = pgl_round(lon_max - lon_min);
1395 /* determine east-west direction */
1396 if (dlon >= 240) {
1397 /* assume that 180th meridian is crossed and swap min/max longitude */
1398 double swap = lon_min; lon_min = lon_max; lon_max = swap;
1399 } else if (dlon > 120) {
1400 /* unclear orientation since delta longitude > 120 */
1401 ereport(ERROR, (
1402 errcode(ERRCODE_DATA_EXCEPTION),
1403 errmsg("can not determine east/west orientation for ebox")
1404 ));
1405 }
1406 /* use boundaries to setup box (and perform checks) */
1407 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1408 /* return result */
1409 PG_RETURN_POINTER(box);
1410 }
1412 /* parse box ("ebox" in SQL) */
1413 /* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
1414 or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
1415 PG_FUNCTION_INFO_V1(pgl_ebox_in);
1416 Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
1417 char *str = PG_GETARG_CSTRING(0); /* input string */
1418 char *str_lower; /* lower case version of input string */
1419 char *strptr; /* current position within string */
1420 int valid; /* number of valid chars */
1421 int done; /* specifies if latitude or longitude was read */
1422 double val; /* temporary variable */
1423 int lat_count = 0; /* count of latitude values parsed */
1424 int lon_count = 0; /* count of longitufde values parsed */
1425 double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */
1426 pgl_box *box; /* return value (to be palloc'ed) */
1427 /* lowercase input */
1428 str_lower = psprintf("%s", str);
1429 for (strptr=str_lower; *strptr; strptr++) {
1430 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1431 }
1432 /* reset reading position to start of (lowercase) string */
1433 strptr = str_lower;
1434 /* check if empty box */
1435 valid = 0;
1436 sscanf(strptr, " empty %n", &valid);
1437 if (valid && strptr[valid] == 0) {
1438 /* allocate and return empty box */
1439 box = (pgl_box *)palloc(sizeof(pgl_box));
1440 pgl_box_set_empty(box);
1441 PG_RETURN_POINTER(box);
1442 }
1443 /* demand four blocks separated by whitespace */
1444 valid = 0;
1445 sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
1446 /* if four blocks separated by whitespace exist, parse those blocks */
1447 if (strptr[valid] == 0) while (strptr[0]) {
1448 /* parse either latitude or longitude (whichever found in input string) */
1449 done = pgl_scan(&strptr, &val, &val);
1450 /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
1451 if (done == PGL_SCAN_LAT) {
1452 if (!lat_count) lat_min = val; else lat_max = val;
1453 lat_count++;
1454 } else if (done == PGL_SCAN_LON) {
1455 if (!lon_count) lon_min = val; else lon_max = val;
1456 lon_count++;
1457 } else {
1458 break;
1459 }
1460 }
1461 /* require end of string, and two latitude and two longitude values */
1462 if (strptr[0] || lat_count != 2 || lon_count != 2) {
1463 ereport(ERROR, (
1464 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1465 errmsg("invalid input syntax for type ebox: \"%s\"", str)
1466 ));
1467 }
1468 /* free lower case string */
1469 pfree(str_lower);
1470 /* order boundaries (maximum greater than minimum) */
1471 if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
1472 if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
1473 /* allocate memory for result */
1474 box = (pgl_box *)palloc(sizeof(pgl_box));
1475 /* set boundaries (and perform checks) */
1476 pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
1477 /* return result */
1478 PG_RETURN_POINTER(box);
1479 }
1481 /* set circle to given latitude, longitude, and radius (including checks) */
1482 static void pgl_ecircle_set_latlon_radius(
1483 pgl_circle *circle, double lat, double lon, double radius
1484 ) {
1485 /* set center point (including checks) */
1486 pgl_epoint_set_latlon(&(circle->center), lat, lon);
1487 /* handle non-positive radius */
1488 if (isnan(radius)) {
1489 ereport(ERROR, (
1490 errcode(ERRCODE_DATA_EXCEPTION),
1491 errmsg("invalid radius for ecircle")
1492 ));
1493 }
1494 if (radius == 0) radius = 0; /* avoids -0 */
1495 else if (radius < 0) {
1496 if (isfinite(radius)) {
1497 ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
1498 }
1499 radius = -INFINITY;
1500 }
1501 /* store radius (round-trip safety is ensured by pgl_print_float) */
1502 circle->radius = radius;
1503 }
1505 /* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
1506 PG_FUNCTION_INFO_V1(pgl_create_ecircle);
1507 Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
1508 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1509 pgl_ecircle_set_latlon_radius(
1510 circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
1511 );
1512 PG_RETURN_POINTER(circle);
1513 }
1515 /* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
1516 PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
1517 Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
1518 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1519 double radius = PG_GETARG_FLOAT8(1);
1520 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1521 /* set latitude, longitude, radius (and perform checks) */
1522 pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
1523 /* return result */
1524 PG_RETURN_POINTER(circle);
1525 }
1527 /* parse circle ("ecircle" in SQL) */
1528 /* format: '[NS]<float> [EW]<float> <float>' */
1529 PG_FUNCTION_INFO_V1(pgl_ecircle_in);
1530 Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
1531 char *str = PG_GETARG_CSTRING(0); /* input string */
1532 char *strptr = str; /* current position within string */
1533 double lat, lon, radius; /* parsed values as double precision flaots */
1534 int valid = 0; /* number of valid chars */
1535 int done = 0; /* stores if latitude and/or longitude was read */
1536 pgl_circle *circle; /* return value (to be palloc'ed) */
1537 /* demand three blocks separated by whitespace */
1538 sscanf(strptr, " %*s %*s %*s %n", &valid);
1539 /* if three blocks separated by whitespace exist, parse those blocks */
1540 if (strptr[valid] == 0) {
1541 /* parse latitude and longitude */
1542 done |= pgl_scan(&strptr, &lat, &lon);
1543 done |= pgl_scan(&strptr, &lat, &lon);
1544 /* parse radius (while incrementing strptr by number of bytes parsed) */
1545 valid = 0;
1546 if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
1547 }
1548 /* require end of string and both latitude and longitude being parsed */
1549 if (strptr[0] || done != PGL_SCAN_LATLON) {
1550 ereport(ERROR, (
1551 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1552 errmsg("invalid input syntax for type ecircle: \"%s\"", str)
1553 ));
1554 }
1555 /* allocate memory for result */
1556 circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1557 /* set latitude, longitude, radius (and perform checks) */
1558 pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
1559 /* return result */
1560 PG_RETURN_POINTER(circle);
1561 }
1563 /* parse cluster ("ecluster" in SQL) */
1564 PG_FUNCTION_INFO_V1(pgl_ecluster_in);
1565 Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
1566 int i;
1567 char *str = PG_GETARG_CSTRING(0); /* input string */
1568 char *str_lower; /* lower case version of input string */
1569 char *strptr; /* pointer to current reading position of input */
1570 int npoints_total = 0; /* total number of points in cluster */
1571 int nentries = 0; /* total number of entries */
1572 pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */
1573 int entries_buflen = 4; /* maximum number of elements in entries array */
1574 int valid; /* number of valid chars processed */
1575 double lat, lon; /* latitude and longitude of parsed point */
1576 int entrytype; /* current entry type */
1577 int npoints; /* number of points in current entry */
1578 pgl_point *points; /* array of pgl_point for pgl_newentry */
1579 int points_buflen; /* maximum number of elements in points array */
1580 int done; /* return value of pgl_scan function */
1581 pgl_cluster *cluster; /* created cluster */
1582 /* lowercase input */
1583 str_lower = psprintf("%s", str);
1584 for (strptr=str_lower; *strptr; strptr++) {
1585 if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
1586 }
1587 /* reset reading position to start of (lowercase) string */
1588 strptr = str_lower;
1589 /* allocate initial buffer for entries */
1590 entries = palloc(entries_buflen * sizeof(pgl_newentry));
1591 /* parse until end of string */
1592 while (strptr[0]) {
1593 /* require previous white-space or closing parenthesis before next token */
1594 if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
1595 goto pgl_ecluster_in_error;
1596 }
1597 /* ignore token "empty" */
1598 valid = 0; sscanf(strptr, " empty %n", &valid);
1599 if (valid) { strptr += valid; continue; }
1600 /* test for "point" token */
1601 valid = 0; sscanf(strptr, " point ( %n", &valid);
1602 if (valid) {
1603 strptr += valid;
1604 entrytype = PGL_ENTRY_POINT;
1605 goto pgl_ecluster_in_type_ok;
1606 }
1607 /* test for "path" token */
1608 valid = 0; sscanf(strptr, " path ( %n", &valid);
1609 if (valid) {
1610 strptr += valid;
1611 entrytype = PGL_ENTRY_PATH;
1612 goto pgl_ecluster_in_type_ok;
1613 }
1614 /* test for "outline" token */
1615 valid = 0; sscanf(strptr, " outline ( %n", &valid);
1616 if (valid) {
1617 strptr += valid;
1618 entrytype = PGL_ENTRY_OUTLINE;
1619 goto pgl_ecluster_in_type_ok;
1620 }
1621 /* test for "polygon" token */
1622 valid = 0; sscanf(strptr, " polygon ( %n", &valid);
1623 if (valid) {
1624 strptr += valid;
1625 entrytype = PGL_ENTRY_POLYGON;
1626 goto pgl_ecluster_in_type_ok;
1627 }
1628 /* error if no valid token found */
1629 goto pgl_ecluster_in_error;
1630 pgl_ecluster_in_type_ok:
1631 /* check if pgl_newentry array needs to grow */
1632 if (nentries == entries_buflen) {
1633 pgl_newentry *newbuf;
1634 entries_buflen *= 2;
1635 newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
1636 memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
1637 pfree(entries);
1638 entries = newbuf;
1639 }
1640 /* reset number of points for current entry */
1641 npoints = 0;
1642 /* allocate array for points */
1643 points_buflen = 4;
1644 points = palloc(points_buflen * sizeof(pgl_point));
1645 /* parse until closing parenthesis */
1646 while (strptr[0] != ')') {
1647 /* error on unexpected end of string */
1648 if (strptr[0] == 0) goto pgl_ecluster_in_error;
1649 /* mark neither latitude nor longitude as read */
1650 done = PGL_SCAN_NONE;
1651 /* require white-space before second, third, etc. point */
1652 if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1653 /* scan latitude (or longitude) */
1654 done |= pgl_scan(&strptr, &lat, &lon);
1655 /* require white-space before second coordinate */
1656 if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
1657 /* scan longitude (or latitude) */
1658 done |= pgl_scan(&strptr, &lat, &lon);
1659 /* error unless both latitude and longitude were parsed */
1660 if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
1661 /* throw error if number of points is too high */
1662 if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
1663 ereport(ERROR, (
1664 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1665 errmsg(
1666 "too many points for ecluster entry (maximum %i)",
1667 PGL_CLUSTER_MAXPOINTS
1668 )
1669 ));
1670 }
1671 /* check if pgl_point array needs to grow */
1672 if (npoints == points_buflen) {
1673 pgl_point *newbuf;
1674 points_buflen *= 2;
1675 newbuf = palloc(points_buflen * sizeof(pgl_point));
1676 memcpy(newbuf, points, npoints * sizeof(pgl_point));
1677 pfree(points);
1678 points = newbuf;
1679 }
1680 /* append point to pgl_point array (includes checks) */
1681 pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
1682 /* increase total number of points */
1683 npoints_total++;
1684 }
1685 /* error if entry has no points */
1686 if (!npoints) goto pgl_ecluster_in_error;
1687 /* entries with one point are automatically of type "point" */
1688 if (npoints == 1) entrytype = PGL_ENTRY_POINT;
1689 /* if entries have more than one point */
1690 else {
1691 /* throw error if entry type is "point" */
1692 if (entrytype == PGL_ENTRY_POINT) {
1693 ereport(ERROR, (
1694 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1695 errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
1696 ));
1697 }
1698 /* coerce outlines and polygons with more than 2 points to be a path */
1699 if (npoints == 2) entrytype = PGL_ENTRY_PATH;
1700 }
1701 /* append entry to pgl_newentry array */
1702 entries[nentries].entrytype = entrytype;
1703 entries[nentries].npoints = npoints;
1704 entries[nentries].points = points;
1705 nentries++;
1706 /* consume closing parenthesis */
1707 strptr++;
1708 /* consume white-space */
1709 while (isspace(strptr[0])) strptr++;
1710 }
1711 /* free lower case string */
1712 pfree(str_lower);
1713 /* create cluster from pgl_newentry array */
1714 cluster = pgl_new_cluster(nentries, entries);
1715 /* free pgl_newentry array */
1716 for (i=0; i<nentries; i++) pfree(entries[i].points);
1717 pfree(entries);
1718 /* set bounding circle of cluster and check east/west orientation */
1719 if (!pgl_finalize_cluster(cluster)) {
1720 ereport(ERROR, (
1721 errcode(ERRCODE_DATA_EXCEPTION),
1722 errmsg("can not determine east/west orientation for ecluster"),
1723 errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
1724 ));
1725 }
1726 /* return cluster */
1727 PG_RETURN_POINTER(cluster);
1728 /* code to throw error */
1729 pgl_ecluster_in_error:
1730 ereport(ERROR, (
1731 errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1732 errmsg("invalid input syntax for type ecluster: \"%s\"", str)
1733 ));
1734 }
1736 /* convert point ("epoint") to string representation */
1737 PG_FUNCTION_INFO_V1(pgl_epoint_out);
1738 Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
1739 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1740 char latstr[PGL_NUMBUFLEN];
1741 char lonstr[PGL_NUMBUFLEN];
1742 pgl_print_lat(latstr, point->lat);
1743 pgl_print_lon(lonstr, point->lon);
1744 PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
1745 }
1747 /* convert box ("ebox") to string representation */
1748 PG_FUNCTION_INFO_V1(pgl_ebox_out);
1749 Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
1750 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
1751 double lon_min = box->lon_min;
1752 double lon_max = box->lon_max;
1753 char lat_min_str[PGL_NUMBUFLEN];
1754 char lat_max_str[PGL_NUMBUFLEN];
1755 char lon_min_str[PGL_NUMBUFLEN];
1756 char lon_max_str[PGL_NUMBUFLEN];
1757 /* return string "empty" if box is set to be empty */
1758 if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
1759 /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
1760 /* (required since pgl_box_in orders the longitude boundaries) */
1761 if (lon_min > lon_max) {
1762 if (lon_min + lon_max >= 0) lon_min -= 360;
1763 else lon_max += 360;
1764 }
1765 /* format and return result */
1766 pgl_print_lat(lat_min_str, box->lat_min);
1767 pgl_print_lat(lat_max_str, box->lat_max);
1768 pgl_print_lon(lon_min_str, lon_min);
1769 pgl_print_lon(lon_max_str, lon_max);
1770 PG_RETURN_CSTRING(psprintf(
1771 "%s %s %s %s",
1772 lat_min_str, lon_min_str, lat_max_str, lon_max_str
1773 ));
1774 }
1776 /* convert circle ("ecircle") to string representation */
1777 PG_FUNCTION_INFO_V1(pgl_ecircle_out);
1778 Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
1779 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
1780 char latstr[PGL_NUMBUFLEN];
1781 char lonstr[PGL_NUMBUFLEN];
1782 char radstr[PGL_NUMBUFLEN];
1783 pgl_print_lat(latstr, circle->center.lat);
1784 pgl_print_lon(lonstr, circle->center.lon);
1785 pgl_print_float(radstr, circle->radius);
1786 PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
1787 }
1789 /* convert cluster ("ecluster") to string representation */
1790 PG_FUNCTION_INFO_V1(pgl_ecluster_out);
1791 Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
1792 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1793 char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */
1794 char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */
1795 char ***strings; /* array of array of strings */
1796 char *string; /* string of current token */
1797 char *res, *resptr; /* result and pointer to current write position */
1798 size_t reslen = 1; /* length of result (init with 1 for terminator) */
1799 int npoints; /* number of points of current entry */
1800 int i, j; /* i: entry, j: point in entry */
1801 /* handle empty clusters */
1802 if (cluster->nentries == 0) {
1803 /* free detoasted cluster (if copy) */
1804 PG_FREE_IF_COPY(cluster, 0);
1805 /* return static result */
1806 PG_RETURN_CSTRING("empty");
1807 }
1808 /* allocate array of array of strings */
1809 strings = palloc(cluster->nentries * sizeof(char **));
1810 /* iterate over all entries in cluster */
1811 for (i=0; i<cluster->nentries; i++) {
1812 /* get number of points in entry */
1813 npoints = cluster->entries[i].npoints;
1814 /* allocate array of strings (one string for each point plus two extra) */
1815 strings[i] = palloc((2 + npoints) * sizeof(char *));
1816 /* determine opening string */
1817 switch (cluster->entries[i].entrytype) {
1818 case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break;
1819 case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break;
1820 case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
1821 case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
1822 default: string = (i==0)?"unknown" :" unknown";
1823 }
1824 /* use opening string as first string in array */
1825 strings[i][0] = string;
1826 /* update result length (for allocating result string later) */
1827 reslen += strlen(string);
1828 /* iterate over all points */
1829 for (j=0; j<npoints; j++) {
1830 /* create string representation of point */
1831 pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
1832 pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
1833 string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
1834 /* copy string pointer to string array */
1835 strings[i][j+1] = string;
1836 /* update result length (for allocating result string later) */
1837 reslen += strlen(string);
1838 }
1839 /* use closing parenthesis as last string in array */
1840 strings[i][npoints+1] = ")";
1841 /* update result length (for allocating result string later) */
1842 reslen++;
1843 }
1844 /* allocate result string */
1845 res = palloc(reslen);
1846 /* set write pointer to begin of result string */
1847 resptr = res;
1848 /* copy strings into result string */
1849 for (i=0; i<cluster->nentries; i++) {
1850 npoints = cluster->entries[i].npoints;
1851 for (j=0; j<npoints+2; j++) {
1852 string = strings[i][j];
1853 strcpy(resptr, string);
1854 resptr += strlen(string);
1855 /* free strings allocated by psprintf */
1856 if (j != 0 && j != npoints+1) pfree(string);
1857 }
1858 /* free array of strings */
1859 pfree(strings[i]);
1860 }
1861 /* free array of array of strings */
1862 pfree(strings);
1863 /* free detoasted cluster (if copy) */
1864 PG_FREE_IF_COPY(cluster, 0);
1865 /* return result */
1866 PG_RETURN_CSTRING(res);
1867 }
1869 /* binary input function for point ("epoint") */
1870 PG_FUNCTION_INFO_V1(pgl_epoint_recv);
1871 Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
1872 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
1873 pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
1874 point->lat = pq_getmsgfloat8(buf);
1875 point->lon = pq_getmsgfloat8(buf);
1876 PG_RETURN_POINTER(point);
1877 }
1879 /* binary input function for box ("ebox") */
1880 PG_FUNCTION_INFO_V1(pgl_ebox_recv);
1881 Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
1882 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
1883 pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
1884 box->lat_min = pq_getmsgfloat8(buf);
1885 box->lat_max = pq_getmsgfloat8(buf);
1886 box->lon_min = pq_getmsgfloat8(buf);
1887 box->lon_max = pq_getmsgfloat8(buf);
1888 PG_RETURN_POINTER(box);
1889 }
1891 /* binary input function for circle ("ecircle") */
1892 PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
1893 Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
1894 StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
1895 pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
1896 circle->center.lat = pq_getmsgfloat8(buf);
1897 circle->center.lon = pq_getmsgfloat8(buf);
1898 circle->radius = pq_getmsgfloat8(buf);
1899 PG_RETURN_POINTER(circle);
1900 }
1902 /* TODO: binary receive function for cluster */
1904 /* binary output function for point ("epoint") */
1905 PG_FUNCTION_INFO_V1(pgl_epoint_send);
1906 Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
1907 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1908 StringInfoData buf;
1909 pq_begintypsend(&buf);
1910 pq_sendfloat8(&buf, point->lat);
1911 pq_sendfloat8(&buf, point->lon);
1912 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1913 }
1915 /* binary output function for box ("ebox") */
1916 PG_FUNCTION_INFO_V1(pgl_ebox_send);
1917 Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
1918 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
1919 StringInfoData buf;
1920 pq_begintypsend(&buf);
1921 pq_sendfloat8(&buf, box->lat_min);
1922 pq_sendfloat8(&buf, box->lat_max);
1923 pq_sendfloat8(&buf, box->lon_min);
1924 pq_sendfloat8(&buf, box->lon_max);
1925 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1926 }
1928 /* binary output function for circle ("ecircle") */
1929 PG_FUNCTION_INFO_V1(pgl_ecircle_send);
1930 Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
1931 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
1932 StringInfoData buf;
1933 pq_begintypsend(&buf);
1934 pq_sendfloat8(&buf, circle->center.lat);
1935 pq_sendfloat8(&buf, circle->center.lon);
1936 pq_sendfloat8(&buf, circle->radius);
1937 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1938 }
1940 /* TODO: binary send functions for cluster */
1942 /* cast point ("epoint") to box ("ebox") */
1943 PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
1944 Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
1945 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1946 pgl_box *box = palloc(sizeof(pgl_box));
1947 box->lat_min = point->lat;
1948 box->lat_max = point->lat;
1949 box->lon_min = point->lon;
1950 box->lon_max = point->lon;
1951 PG_RETURN_POINTER(box);
1952 }
1954 /* cast point ("epoint") to circle ("ecircle") */
1955 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
1956 Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
1957 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1958 pgl_circle *circle = palloc(sizeof(pgl_box));
1959 circle->center = *point;
1960 circle->radius = 0;
1961 PG_RETURN_POINTER(circle);
1962 }
1964 /* cast point ("epoint") to cluster ("ecluster") */
1965 PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
1966 Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
1967 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
1968 pgl_newentry entry;
1969 entry.entrytype = PGL_ENTRY_POINT;
1970 entry.npoints = 1;
1971 entry.points = point;
1972 PG_RETURN_POINTER(pgl_new_cluster(1, &entry));
1973 }
1975 /* cast box ("ebox") to cluster ("ecluster") */
1976 #define pgl_ebox_to_ecluster_macro(i, a, b) \
1977 entries[i].entrytype = PGL_ENTRY_POLYGON; \
1978 entries[i].npoints = 4; \
1979 entries[i].points = points[i]; \
1980 points[i][0].lat = box->lat_min; \
1981 points[i][0].lon = (a); \
1982 points[i][1].lat = box->lat_min; \
1983 points[i][1].lon = (b); \
1984 points[i][2].lat = box->lat_max; \
1985 points[i][2].lon = (b); \
1986 points[i][3].lat = box->lat_max; \
1987 points[i][3].lon = (a);
1988 PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
1989 Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
1990 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
1991 double lon, dlon;
1992 int nentries;
1993 pgl_newentry entries[3];
1994 pgl_point points[3][4];
1995 if (box->lat_min > box->lat_max) {
1996 nentries = 0;
1997 } else if (box->lon_min > box->lon_max) {
1998 if (box->lon_min < 0) {
1999 lon = pgl_round((box->lon_min + 180) / 2.0);
2000 nentries = 3;
2001 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2002 pgl_ebox_to_ecluster_macro(1, lon, 180);
2003 pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
2004 } else if (box->lon_max > 0) {
2005 lon = pgl_round((box->lon_max - 180) / 2.0);
2006 nentries = 3;
2007 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2008 pgl_ebox_to_ecluster_macro(1, -180, lon);
2009 pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
2010 } else {
2011 nentries = 2;
2012 pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
2013 pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
2014 }
2015 } else {
2016 dlon = pgl_round(box->lon_max - box->lon_min);
2017 if (dlon < 180) {
2018 nentries = 1;
2019 pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
2020 } else {
2021 lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
2022 if (
2023 pgl_round(lon - box->lon_min) < 180 &&
2024 pgl_round(box->lon_max - lon) < 180
2025 ) {
2026 nentries = 2;
2027 pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
2028 pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
2029 } else {
2030 nentries = 3;
2031 pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
2032 pgl_ebox_to_ecluster_macro(1, -60, 60);
2033 pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
2034 }
2035 }
2036 }
2037 PG_RETURN_POINTER(pgl_new_cluster(nentries, entries));
2038 }
2040 /* extract latitude from point ("epoint") */
2041 PG_FUNCTION_INFO_V1(pgl_epoint_lat);
2042 Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
2043 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
2044 }
2046 /* extract longitude from point ("epoint") */
2047 PG_FUNCTION_INFO_V1(pgl_epoint_lon);
2048 Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
2049 PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
2050 }
2052 /* extract minimum latitude from box ("ebox") */
2053 PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
2054 Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
2055 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
2056 }
2058 /* extract maximum latitude from box ("ebox") */
2059 PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
2060 Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
2061 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
2062 }
2064 /* extract minimum longitude from box ("ebox") */
2065 PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
2066 Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
2067 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
2068 }
2070 /* extract maximum longitude from box ("ebox") */
2071 PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
2072 Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
2073 PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
2074 }
2076 /* extract center point from circle ("ecircle") */
2077 PG_FUNCTION_INFO_V1(pgl_ecircle_center);
2078 Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
2079 PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
2080 }
2082 /* extract radius from circle ("ecircle") */
2083 PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
2084 Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
2085 PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
2086 }
2088 /* check if point is inside box (overlap operator "&&") in SQL */
2089 PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
2090 Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
2091 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2092 pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
2093 PG_RETURN_BOOL(pgl_point_in_box(point, box));
2094 }
2096 /* check if point is inside circle (overlap operator "&&") in SQL */
2097 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
2098 Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
2099 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2100 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2101 PG_RETURN_BOOL(
2102 pgl_distance(
2103 point->lat, point->lon,
2104 circle->center.lat, circle->center.lon
2105 ) <= circle->radius
2106 );
2107 }
2109 /* check if point is inside cluster (overlap operator "&&") in SQL */
2110 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
2111 Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
2112 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2113 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2114 bool retval = pgl_point_in_cluster(point, cluster);
2115 PG_FREE_IF_COPY(cluster, 1);
2116 PG_RETURN_BOOL(retval);
2117 }
2119 /* check if two boxes overlap (overlap operator "&&") in SQL */
2120 PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
2121 Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
2122 pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
2123 pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
2124 PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
2125 }
2127 /* check if two circles overlap (overlap operator "&&") in SQL */
2128 PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
2129 Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
2130 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2131 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2132 PG_RETURN_BOOL(
2133 pgl_distance(
2134 circle1->center.lat, circle1->center.lon,
2135 circle2->center.lat, circle2->center.lon
2136 ) <= circle1->radius + circle2->radius
2137 );
2138 }
2140 /* check if circle and cluster overlap (overlap operator "&&") in SQL */
2141 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
2142 Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
2143 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2144 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2145 bool retval = (
2146 pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
2147 );
2148 PG_FREE_IF_COPY(cluster, 1);
2149 PG_RETURN_BOOL(retval);
2150 }
2152 /* calculate distance between two points ("<->" operator) in SQL */
2153 PG_FUNCTION_INFO_V1(pgl_epoint_distance);
2154 Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
2155 pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
2156 pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
2157 PG_RETURN_FLOAT8(pgl_distance(
2158 point1->lat, point1->lon, point2->lat, point2->lon
2159 ));
2160 }
2162 /* calculate point to circle distance ("<->" operator) in SQL */
2163 PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
2164 Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
2165 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2166 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
2167 double distance = pgl_distance(
2168 point->lat, point->lon, circle->center.lat, circle->center.lon
2169 ) - circle->radius;
2170 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2171 }
2173 /* calculate point to cluster distance ("<->" operator) in SQL */
2174 PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
2175 Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
2176 pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
2177 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2178 double distance = pgl_point_cluster_distance(point, cluster);
2179 PG_FREE_IF_COPY(cluster, 1);
2180 PG_RETURN_FLOAT8(distance);
2181 }
2183 /* calculate distance between two circles ("<->" operator) in SQL */
2184 PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
2185 Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
2186 pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
2187 pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
2188 double distance = pgl_distance(
2189 circle1->center.lat, circle1->center.lon,
2190 circle2->center.lat, circle2->center.lon
2191 ) - (circle1->radius + circle2->radius);
2192 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2193 }
2195 /* calculate circle to cluster distance ("<->" operator) in SQL */
2196 PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
2197 Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
2198 pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
2199 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2200 double distance = (
2201 pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
2202 );
2203 PG_FREE_IF_COPY(cluster, 1);
2204 PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
2205 }
2208 /*-----------------------------------------------------------*
2209 * B-tree comparison operators and index support functions *
2210 *-----------------------------------------------------------*/
2212 /* macro for a B-tree operator (without detoasting) */
2213 #define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
2214 PG_FUNCTION_INFO_V1(func); \
2215 Datum func(PG_FUNCTION_ARGS) { \
2216 type *a = (type *)PG_GETARG_POINTER(0); \
2217 type *b = (type *)PG_GETARG_POINTER(1); \
2218 PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
2219 }
2221 /* macro for a B-tree comparison function (without detoasting) */
2222 #define PGL_BTREE_CMP(func, type, cmpfunc) \
2223 PG_FUNCTION_INFO_V1(func); \
2224 Datum func(PG_FUNCTION_ARGS) { \
2225 type *a = (type *)PG_GETARG_POINTER(0); \
2226 type *b = (type *)PG_GETARG_POINTER(1); \
2227 PG_RETURN_INT32(cmpfunc(a, b)); \
2228 }
2230 /* macro for a B-tree operator (with detoasting) */
2231 #define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
2232 PG_FUNCTION_INFO_V1(func); \
2233 Datum func(PG_FUNCTION_ARGS) { \
2234 bool res; \
2235 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2236 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2237 res = cmpfunc(a, b) oper 0; \
2238 PG_FREE_IF_COPY(a, 0); \
2239 PG_FREE_IF_COPY(b, 1); \
2240 PG_RETURN_BOOL(res); \
2241 }
2243 /* macro for a B-tree comparison function (with detoasting) */
2244 #define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
2245 PG_FUNCTION_INFO_V1(func); \
2246 Datum func(PG_FUNCTION_ARGS) { \
2247 int32_t res; \
2248 type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
2249 type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
2250 res = cmpfunc(a, b); \
2251 PG_FREE_IF_COPY(a, 0); \
2252 PG_FREE_IF_COPY(b, 1); \
2253 PG_RETURN_INT32(res); \
2254 }
2256 /* B-tree operators and comparison function for point */
2257 PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
2258 PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
2259 PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
2260 PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
2261 PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
2262 PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
2263 PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
2265 /* B-tree operators and comparison function for box */
2266 PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
2267 PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
2268 PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
2269 PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
2270 PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
2271 PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
2272 PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
2274 /* B-tree operators and comparison function for circle */
2275 PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
2276 PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
2277 PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
2278 PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
2279 PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
2280 PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
2281 PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
2284 /*--------------------------------*
2285 * GiST index support functions *
2286 *--------------------------------*/
2288 /* GiST "consistent" support function */
2289 PG_FUNCTION_INFO_V1(pgl_gist_consistent);
2290 Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
2291 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2292 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
2293 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
2294 bool *recheck = (bool *)PG_GETARG_POINTER(4);
2295 /* demand recheck because index and query methods are lossy */
2296 *recheck = true;
2297 /* strategy number 11: equality of two points */
2298 if (strategy == 11) {
2299 /* query datum is another point */
2300 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2301 /* convert other point to key */
2302 pgl_pointkey querykey;
2303 pgl_point_to_key(query, querykey);
2304 /* return true if both keys overlap */
2305 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2306 }
2307 /* strategy number 13: equality of two circles */
2308 if (strategy == 13) {
2309 /* query datum is another circle */
2310 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2311 /* convert other circle to key */
2312 pgl_areakey querykey;
2313 pgl_circle_to_key(query, querykey);
2314 /* return true if both keys overlap */
2315 PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
2316 }
2317 /* for all remaining strategies, keys on empty objects produce no match */
2318 /* (check necessary because query radius may be infinite) */
2319 if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
2320 /* strategy number 21: overlapping with point */
2321 if (strategy == 21) {
2322 /* query datum is a point */
2323 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2324 /* return true if estimated distance (allowed to be smaller than real
2325 distance) between index key and point is zero */
2326 PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
2327 }
2328 /* strategy number 22: (point) overlapping with box */
2329 if (strategy == 22) {
2330 /* query datum is a box */
2331 pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
2332 /* determine bounding box of indexed key */
2333 pgl_box keybox;
2334 pgl_key_to_box(key, &keybox);
2335 /* return true if query box overlaps with bounding box of indexed key */
2336 PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
2337 }
2338 /* strategy number 23: overlapping with circle */
2339 if (strategy == 23) {
2340 /* query datum is a circle */
2341 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2342 /* return true if estimated distance (allowed to be smaller than real
2343 distance) between index key and circle center is smaller than radius */
2344 PG_RETURN_BOOL(
2345 pgl_estimate_key_distance(key, &(query->center)) <= query->radius
2346 );
2347 }
2348 /* strategy number 24: overlapping with cluster */
2349 if (strategy == 24) {
2350 bool retval; /* return value */
2351 /* query datum is a cluster */
2352 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2353 /* return true if estimated distance (allowed to be smaller than real
2354 distance) between index key and circle center is smaller than radius */
2355 retval = (
2356 pgl_estimate_key_distance(key, &(query->bounding.center)) <=
2357 query->bounding.radius
2358 );
2359 PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */
2360 PG_RETURN_BOOL(retval);
2361 }
2362 /* throw error for any unknown strategy number */
2363 elog(ERROR, "unrecognized strategy number: %d", strategy);
2364 }
2366 /* GiST "union" support function */
2367 PG_FUNCTION_INFO_V1(pgl_gist_union);
2368 Datum pgl_gist_union(PG_FUNCTION_ARGS) {
2369 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2370 pgl_keyptr out; /* return value (to be palloc'ed) */
2371 int i;
2372 /* determine key size */
2373 size_t keysize = PGL_KEY_IS_AREAKEY(
2374 (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
2375 ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
2376 /* begin with first key as result */
2377 out = palloc(keysize);
2378 memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
2379 /* unite current result with second, third, etc. key */
2380 for (i=1; i<entryvec->n; i++) {
2381 pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
2382 }
2383 /* return result */
2384 PG_RETURN_POINTER(out);
2385 }
2387 /* GiST "compress" support function for indicis on points */
2388 PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
2389 Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
2390 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2391 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2392 /* only transform new leaves */
2393 if (entry->leafkey) {
2394 /* get point to be transformed */
2395 pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
2396 /* allocate memory for key */
2397 pgl_keyptr key = palloc(sizeof(pgl_pointkey));
2398 /* transform point to key */
2399 pgl_point_to_key(point, key);
2400 /* create new GISTENTRY structure as return value */
2401 retval = palloc(sizeof(GISTENTRY));
2402 gistentryinit(
2403 *retval, PointerGetDatum(key),
2404 entry->rel, entry->page, entry->offset, FALSE
2405 );
2406 } else {
2407 /* inner nodes have already been transformed */
2408 retval = entry;
2409 }
2410 /* return pointer to old or new GISTENTRY structure */
2411 PG_RETURN_POINTER(retval);
2412 }
2414 /* GiST "compress" support function for indicis on circles */
2415 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
2416 Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
2417 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2418 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2419 /* only transform new leaves */
2420 if (entry->leafkey) {
2421 /* get circle to be transformed */
2422 pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
2423 /* allocate memory for key */
2424 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2425 /* transform circle to key */
2426 pgl_circle_to_key(circle, key);
2427 /* create new GISTENTRY structure as return value */
2428 retval = palloc(sizeof(GISTENTRY));
2429 gistentryinit(
2430 *retval, PointerGetDatum(key),
2431 entry->rel, entry->page, entry->offset, FALSE
2432 );
2433 } else {
2434 /* inner nodes have already been transformed */
2435 retval = entry;
2436 }
2437 /* return pointer to old or new GISTENTRY structure */
2438 PG_RETURN_POINTER(retval);
2439 }
2441 /* GiST "compress" support function for indices on clusters */
2442 PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
2443 Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
2444 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
2445 GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
2446 /* only transform new leaves */
2447 if (entry->leafkey) {
2448 /* get cluster to be transformed (detoasting necessary!) */
2449 pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
2450 /* allocate memory for key */
2451 pgl_keyptr key = palloc(sizeof(pgl_areakey));
2452 /* transform cluster to key */
2453 pgl_circle_to_key(&(cluster->bounding), key);
2454 /* create new GISTENTRY structure as return value */
2455 retval = palloc(sizeof(GISTENTRY));
2456 gistentryinit(
2457 *retval, PointerGetDatum(key),
2458 entry->rel, entry->page, entry->offset, FALSE
2459 );
2460 /* free detoasted datum */
2461 if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
2462 } else {
2463 /* inner nodes have already been transformed */
2464 retval = entry;
2465 }
2466 /* return pointer to old or new GISTENTRY structure */
2467 PG_RETURN_POINTER(retval);
2468 }
2470 /* GiST "decompress" support function for indices */
2471 PG_FUNCTION_INFO_V1(pgl_gist_decompress);
2472 Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
2473 /* return passed pointer without transformation */
2474 PG_RETURN_POINTER(PG_GETARG_POINTER(0));
2475 }
2477 /* GiST "penalty" support function */
2478 PG_FUNCTION_INFO_V1(pgl_gist_penalty);
2479 Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
2480 GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
2481 GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
2482 float *penalty = (float *)PG_GETARG_POINTER(2);
2483 /* get original key and key to insert */
2484 pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
2485 pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
2486 /* copy original key */
2487 union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
2488 if (PGL_KEY_IS_AREAKEY(orig)) {
2489 memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
2490 } else {
2491 memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
2492 }
2493 /* calculate union of both keys */
2494 pgl_unite_keys((pgl_keyptr)&union_key, new);
2495 /* penalty equal to reduction of key length (logarithm of added area) */
2496 /* (return value by setting referenced value and returning pointer) */
2497 *penalty = (
2498 PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
2499 );
2500 PG_RETURN_POINTER(penalty);
2501 }
2503 /* GiST "picksplit" support function */
2504 PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
2505 Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
2506 GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
2507 GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
2508 OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */
2509 union {
2510 pgl_pointkey pointkey;
2511 pgl_areakey areakey;
2512 } union_all; /* union of all keys (to be calculated from scratch)
2513 (later cut in half) */
2514 int is_areakey = PGL_KEY_IS_AREAKEY(
2515 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
2516 );
2517 int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
2518 pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */
2519 pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */
2520 pgl_keyptr key; /* current key to be processed */
2521 /* allocate memory for array of left and right keys, set counts to zero */
2522 v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2523 v->spl_nleft = 0;
2524 v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
2525 v->spl_nright = 0;
2526 /* calculate union of all keys from scratch */
2527 memcpy(
2528 (pgl_keyptr)&union_all,
2529 (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
2530 keysize
2531 );
2532 for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
2533 pgl_unite_keys(
2534 (pgl_keyptr)&union_all,
2535 (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
2536 );
2537 }
2538 /* check if trivial split is necessary due to exhausted key length */
2539 /* (Note: keys for empty objects must have node depth set to maximum) */
2540 if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
2541 is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
2542 )) {
2543 /* half of all keys go left */
2544 for (
2545 i=FirstOffsetNumber;
2546 i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2547 i=OffsetNumberNext(i)
2548 ) {
2549 /* pointer to current key */
2550 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2551 /* update unionL */
2552 /* check if key is first key that goes left */
2553 if (!v->spl_nleft) {
2554 /* first key that goes left is just copied to unionL */
2555 memcpy(unionL, key, keysize);
2556 } else {
2557 /* unite current value and next key */
2558 pgl_unite_keys(unionL, key);
2559 }
2560 /* append offset number to list of keys that go left */
2561 v->spl_left[v->spl_nleft++] = i;
2562 }
2563 /* other half goes right */
2564 for (
2565 i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
2566 i<entryvec->n;
2567 i=OffsetNumberNext(i)
2568 ) {
2569 /* pointer to current key */
2570 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2571 /* update unionR */
2572 /* check if key is first key that goes right */
2573 if (!v->spl_nright) {
2574 /* first key that goes right is just copied to unionR */
2575 memcpy(unionR, key, keysize);
2576 } else {
2577 /* unite current value and next key */
2578 pgl_unite_keys(unionR, key);
2579 }
2580 /* append offset number to list of keys that go right */
2581 v->spl_right[v->spl_nright++] = i;
2582 }
2583 }
2584 /* otherwise, a non-trivial split is possible */
2585 else {
2586 /* cut covered area in half */
2587 /* (union_all then refers to area of keys that go left) */
2588 /* check if union of all keys covers empty and non-empty objects */
2589 if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
2590 /* if yes, split into empty and non-empty objects */
2591 pgl_key_set_empty((pgl_keyptr)&union_all);
2592 } else {
2593 /* otherwise split by next bit */
2594 ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
2595 /* NOTE: type bit conserved */
2596 }
2597 /* determine for each key if it goes left or right */
2598 for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
2599 /* pointer to current key */
2600 key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
2601 /* keys within one half of the area go left */
2602 if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
2603 /* update unionL */
2604 /* check if key is first key that goes left */
2605 if (!v->spl_nleft) {
2606 /* first key that goes left is just copied to unionL */
2607 memcpy(unionL, key, keysize);
2608 } else {
2609 /* unite current value of unionL and processed key */
2610 pgl_unite_keys(unionL, key);
2611 }
2612 /* append offset number to list of keys that go left */
2613 v->spl_left[v->spl_nleft++] = i;
2614 }
2615 /* the other keys go right */
2616 else {
2617 /* update unionR */
2618 /* check if key is first key that goes right */
2619 if (!v->spl_nright) {
2620 /* first key that goes right is just copied to unionR */
2621 memcpy(unionR, key, keysize);
2622 } else {
2623 /* unite current value of unionR and processed key */
2624 pgl_unite_keys(unionR, key);
2625 }
2626 /* append offset number to list of keys that go right */
2627 v->spl_right[v->spl_nright++] = i;
2628 }
2629 }
2630 }
2631 /* store unions in return value */
2632 v->spl_ldatum = PointerGetDatum(unionL);
2633 v->spl_rdatum = PointerGetDatum(unionR);
2634 /* return all results */
2635 PG_RETURN_POINTER(v);
2636 }
2638 /* GiST "same"/"equal" support function */
2639 PG_FUNCTION_INFO_V1(pgl_gist_same);
2640 Datum pgl_gist_same(PG_FUNCTION_ARGS) {
2641 pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
2642 pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
2643 bool *boolptr = (bool *)PG_GETARG_POINTER(2);
2644 /* two keys are equal if they are binary equal */
2645 /* (return result by setting referenced boolean and returning pointer) */
2646 *boolptr = !memcmp(
2647 key1,
2648 key2,
2649 PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
2650 );
2651 PG_RETURN_POINTER(boolptr);
2652 }
2654 /* GiST "distance" support function */
2655 PG_FUNCTION_INFO_V1(pgl_gist_distance);
2656 Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
2657 GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
2658 pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
2659 StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
2660 bool *recheck = (bool *)PG_GETARG_POINTER(4);
2661 double distance; /* return value */
2662 /* demand recheck because distance is just an estimation */
2663 /* (real distance may be bigger) */
2664 *recheck = true;
2665 /* strategy number 31: distance to point */
2666 if (strategy == 31) {
2667 /* query datum is a point */
2668 pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
2669 /* use pgl_estimate_pointkey_distance() function to compute result */
2670 distance = pgl_estimate_key_distance(key, query);
2671 /* avoid infinity (reserved!) */
2672 if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
2673 /* return result */
2674 PG_RETURN_FLOAT8(distance);
2675 }
2676 /* strategy number 33: distance to circle */
2677 if (strategy == 33) {
2678 /* query datum is a circle */
2679 pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
2680 /* estimate distance to circle center and substract circle radius */
2681 distance = (
2682 pgl_estimate_key_distance(key, &(query->center)) - query->radius
2683 );
2684 /* convert non-positive values to zero and avoid infinity (reserved!) */
2685 if (distance <= 0) distance = 0;
2686 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
2687 /* return result */
2688 PG_RETURN_FLOAT8(distance);
2689 }
2690 /* strategy number 34: distance to cluster */
2691 if (strategy == 34) {
2692 /* query datum is a cluster */
2693 pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2694 /* estimate distance to bounding center and substract bounding radius */
2695 distance = (
2696 pgl_estimate_key_distance(key, &(query->bounding.center)) -
2697 query->bounding.radius
2698 );
2699 /* convert non-positive values to zero and avoid infinity (reserved!) */
2700 if (distance <= 0) distance = 0;
2701 else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
2702 /* free detoasted cluster (if copy) */
2703 PG_FREE_IF_COPY(query, 1);
2704 /* return result */
2705 PG_RETURN_FLOAT8(distance);
2706 }
2707 /* throw error for any unknown strategy number */
2708 elog(ERROR, "unrecognized strategy number: %d", strategy);
2709 }