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