pgLatLon
diff latlon-v0010.c @ 69:882b77aee653
Prepared file names for version 0.15
author | jbe |
---|---|
date | Wed Feb 12 11:08:37 2020 +0100 (2020-02-12) |
parents | latlon-v0009.c@a5b8024ef5bc |
children | d06b066fb8ad |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/latlon-v0010.c Wed Feb 12 11:08:37 2020 +0100 1.3 @@ -0,0 +1,3352 @@ 1.4 + 1.5 +/*-------------* 1.6 + * C prelude * 1.7 + *-------------*/ 1.8 + 1.9 +#include "postgres.h" 1.10 +#include "fmgr.h" 1.11 +#include "libpq/pqformat.h" 1.12 +#include "access/gist.h" 1.13 +#include "access/stratnum.h" 1.14 +#include "utils/array.h" 1.15 +#include <limits.h> 1.16 +#include <math.h> 1.17 + 1.18 +#ifdef PG_MODULE_MAGIC 1.19 +PG_MODULE_MAGIC; 1.20 +#endif 1.21 + 1.22 +#if INT_MAX < 2147483647 1.23 +#error Expected int type to be at least 32 bit wide 1.24 +#endif 1.25 + 1.26 + 1.27 +/*---------------------------------* 1.28 + * distance calculation on earth * 1.29 + * (using WGS-84 spheroid) * 1.30 + *---------------------------------*/ 1.31 + 1.32 +/* WGS-84 spheroid with following parameters: 1.33 + semi-major axis a = 6378137 1.34 + semi-minor axis b = a * (1 - 1/298.257223563) 1.35 + estimated diameter = 2 * (2*a+b)/3 1.36 +*/ 1.37 +#define PGL_SPHEROID_A 6378137.0 /* semi major axis */ 1.38 +#define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */ 1.39 +#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F)) 1.40 +#define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \ 1.41 + PGL_SPHEROID_B * PGL_SPHEROID_B ) / \ 1.42 + ( PGL_SPHEROID_A * PGL_SPHEROID_A ) ) 1.43 +#define PGL_SUBEPS2 (1.0-PGL_EPS2) 1.44 +#define PGL_RADIUS ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0) 1.45 +#define PGL_DIAMETER (2.0 * PGL_RADIUS) 1.46 +#define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */ 1.47 +#define PGL_MAXDIST (PGL_RADIUS * M_PI) /* maximum distance */ 1.48 +#define PGL_FADELIMIT (PGL_MAXDIST / 3.0) /* 1/6 circumference */ 1.49 + 1.50 +/* calculate distance between two points on earth (given in degrees) */ 1.51 +static inline double pgl_distance( 1.52 + double lat1, double lon1, double lat2, double lon2 1.53 +) { 1.54 + float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin; 1.55 + float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t; 1.56 + /* normalize delta longitude (lon2 > 0 && lon1 = 0) */ 1.57 + /* lon1 = 0 (not used anymore) */ 1.58 + lon2 = fabs(lon2-lon1); 1.59 + /* convert to radians (first divide, then multiply) */ 1.60 + lat1 = (lat1 / 180.0) * M_PI; 1.61 + lat2 = (lat2 / 180.0) * M_PI; 1.62 + lon2 = (lon2 / 180.0) * M_PI; 1.63 + /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point 1.64 + operations (lon2 >= lon1 is already ensured in a previous step) */ 1.65 + if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; } 1.66 + /* calculate 3d coordinates on scaled ellipsoid which has an average diameter 1.67 + of 1.0 */ 1.68 + lat1cos = cos(lat1); lat1sin = sin(lat1); 1.69 + lat2cos = cos(lat2); lat2sin = sin(lat2); 1.70 + lon2cos = cos(lon2); lon2sin = sin(lon2); 1.71 + nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin); 1.72 + nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin); 1.73 + x1 = nphi1 * lat1cos; 1.74 + z1 = nphi1 * PGL_SUBEPS2 * lat1sin; 1.75 + x2 = nphi2 * lat2cos * lon2cos; 1.76 + y2 = nphi2 * lat2cos * lon2sin; 1.77 + z2 = nphi2 * PGL_SUBEPS2 * lat2sin; 1.78 + /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */ 1.79 + g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1)); 1.80 + /* convert tunnel distance through scaled ellipsoid to approximated surface 1.81 + distance on original ellipsoid */ 1.82 + if (g > 1.0) g = 1.0; 1.83 + s = PGL_DIAMETER * asin(g); 1.84 + /* return result only if small enough to be precise (less than 1/3 of 1.85 + maximum possible distance) */ 1.86 + if (s <= PGL_FADELIMIT) return s; 1.87 + /* calculate tunnel distance to antipodal point through scaled ellipsoid */ 1.88 + g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1)); 1.89 + /* convert tunnel distance to antipodal point through scaled ellipsoid to 1.90 + approximated surface distance to antipodal point on original ellipsoid */ 1.91 + if (g > 1.0) g = 1.0; 1.92 + t = PGL_DIAMETER * asin(g); 1.93 + /* surface distance between original points can now be approximated by 1.94 + substracting antipodal distance from maximum possible distance; 1.95 + return result only if small enough (less than 1/3 of maximum possible 1.96 + distance) */ 1.97 + if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t; 1.98 + /* otherwise crossfade direct and antipodal result to ensure monotonicity */ 1.99 + return ( 1.100 + (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) / 1.101 + (s + t - 2*PGL_FADELIMIT) 1.102 + ); 1.103 +} 1.104 + 1.105 +/* finite distance that can not be reached on earth */ 1.106 +#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST) 1.107 + 1.108 + 1.109 +/*--------------------------------* 1.110 + * simple geographic data types * 1.111 + *--------------------------------*/ 1.112 + 1.113 +/* point on earth given by latitude and longitude in degrees */ 1.114 +/* (type "epoint" in SQL) */ 1.115 +typedef struct { 1.116 + double lat; /* between -90 and 90 (both inclusive) */ 1.117 + double lon; /* between -180 and 180 (both inclusive) */ 1.118 +} pgl_point; 1.119 + 1.120 +/* box delimited by two parallels and two meridians (all in degrees) */ 1.121 +/* (type "ebox" in SQL) */ 1.122 +typedef struct { 1.123 + double lat_min; /* between -90 and 90 (both inclusive) */ 1.124 + double lat_max; /* between -90 and 90 (both inclusive) */ 1.125 + double lon_min; /* between -180 and 180 (both inclusive) */ 1.126 + double lon_max; /* between -180 and 180 (both inclusive) */ 1.127 + /* if lat_min > lat_max, then box is empty */ 1.128 + /* if lon_min > lon_max, then 180th meridian is crossed */ 1.129 +} pgl_box; 1.130 + 1.131 +/* circle on earth surface (for radial searches with fixed radius) */ 1.132 +/* (type "ecircle" in SQL) */ 1.133 +typedef struct { 1.134 + pgl_point center; 1.135 + double radius; /* positive (including +0 but excluding -0), or -INFINITY */ 1.136 + /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point), 1.137 + zero radius (0) denotes a single point, 1.138 + a finite radius (0 < radius < INFINITY) denotes a filled circle, and 1.139 + a radius of INFINITY is valid and means complete coverage of earth. */ 1.140 +} pgl_circle; 1.141 + 1.142 + 1.143 +/*----------------------------------* 1.144 + * geographic "cluster" data type * 1.145 + *----------------------------------*/ 1.146 + 1.147 +/* A cluster is a collection of points, paths, outlines, and polygons. If two 1.148 + polygons in a cluster overlap, the area covered by both polygons does not 1.149 + belong to the cluster. This way, a cluster can be used to describe complex 1.150 + shapes like polygons with holes. Outlines are non-filled polygons. Paths are 1.151 + open by default (i.e. the last point in the list is not connected with the 1.152 + first point in the list). Note that each outline or polygon in a cluster 1.153 + must cover a longitude range of less than 180 degrees to avoid ambiguities. 1.154 + Areas which are larger may be split into multiple polygons. */ 1.155 + 1.156 +/* maximum number of points in a cluster */ 1.157 +/* (limited to avoid integer overflows, e.g. when allocating memory) */ 1.158 +#define PGL_CLUSTER_MAXPOINTS 16777216 1.159 + 1.160 +/* types of cluster entries */ 1.161 +#define PGL_ENTRY_POINT 1 /* a point */ 1.162 +#define PGL_ENTRY_PATH 2 /* a path from first point to last point */ 1.163 +#define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */ 1.164 +#define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */ 1.165 + 1.166 +/* Entries of a cluster are described by two different structs: pgl_newentry 1.167 + and pgl_entry. The first is used only during construction of a cluster, the 1.168 + second is used in all other cases (e.g. when reading clusters from the 1.169 + database, performing operations, etc). */ 1.170 + 1.171 +/* entry for new geographic cluster during construction of that cluster */ 1.172 +typedef struct { 1.173 + int32_t entrytype; 1.174 + int32_t npoints; 1.175 + pgl_point *points; /* pointer to an array of points (pgl_point) */ 1.176 +} pgl_newentry; 1.177 + 1.178 +/* entry of geographic cluster */ 1.179 +typedef struct { 1.180 + int32_t entrytype; /* type of entry: point, path, outline, polygon */ 1.181 + int32_t npoints; /* number of stored points (set to 1 for point entry) */ 1.182 + int32_t offset; /* offset of pgl_point array from cluster base address */ 1.183 + /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */ 1.184 +} pgl_entry; 1.185 + 1.186 +/* geographic cluster which is a collection of points, (open) paths, polygons, 1.187 + and outlines (non-filled polygons) */ 1.188 +typedef struct { 1.189 + char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */ 1.190 + int32_t nentries; /* number of stored points */ 1.191 + pgl_circle bounding; /* bounding circle */ 1.192 + /* Note: bounding circle ensures alignment of pgl_cluster for points */ 1.193 + pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */ 1.194 +} pgl_cluster; 1.195 + 1.196 +/* macro to determine memory alignment of points */ 1.197 +/* (needed to store pgl_point array after entries in pgl_cluster) */ 1.198 +typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment; 1.199 +#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned) 1.200 + 1.201 +/* macro to extract a pointer to the array of points of a cluster entry */ 1.202 +#define PGL_ENTRY_POINTS(cluster, idx) \ 1.203 + ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset)) 1.204 + 1.205 +/* convert pgl_newentry array to pgl_cluster */ 1.206 +/* NOTE: requires pgl_finalize_cluster to be called to finalize result */ 1.207 +static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) { 1.208 + int i; /* index of current entry */ 1.209 + int npoints = 0; /* number of points in whole cluster */ 1.210 + int entry_npoints; /* number of points in current entry */ 1.211 + int points_offset = PGL_POINT_ALIGNMENT * ( 1.212 + ( offsetof(pgl_cluster, entries) + 1.213 + nentries * sizeof(pgl_entry) + 1.214 + PGL_POINT_ALIGNMENT - 1 1.215 + ) / PGL_POINT_ALIGNMENT 1.216 + ); /* offset of pgl_point array from base address (considering alignment) */ 1.217 + pgl_cluster *cluster; /* new cluster to be returned */ 1.218 + /* determine total number of points */ 1.219 + for (i=0; i<nentries; i++) npoints += entries[i].npoints; 1.220 + /* allocate memory for cluster (including entries and points) */ 1.221 + cluster = palloc(points_offset + npoints * sizeof(pgl_point)); 1.222 + /* re-count total number of points to determine offset for each entry */ 1.223 + npoints = 0; 1.224 + /* copy entries and points */ 1.225 + for (i=0; i<nentries; i++) { 1.226 + /* determine number of points in entry */ 1.227 + entry_npoints = entries[i].npoints; 1.228 + /* copy entry */ 1.229 + cluster->entries[i].entrytype = entries[i].entrytype; 1.230 + cluster->entries[i].npoints = entry_npoints; 1.231 + /* calculate offset (in bytes) of pgl_point array */ 1.232 + cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point); 1.233 + /* copy points */ 1.234 + memcpy( 1.235 + PGL_ENTRY_POINTS(cluster, i), 1.236 + entries[i].points, 1.237 + entry_npoints * sizeof(pgl_point) 1.238 + ); 1.239 + /* update total number of points processed */ 1.240 + npoints += entry_npoints; 1.241 + } 1.242 + /* set number of entries in cluster */ 1.243 + cluster->nentries = nentries; 1.244 + /* set PostgreSQL header for variable sized data */ 1.245 + SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point)); 1.246 + /* return newly created cluster */ 1.247 + return cluster; 1.248 +} 1.249 + 1.250 + 1.251 +/*----------------------------------------------* 1.252 + * Geographic point with integer sample count * 1.253 + * (needed for fair distance calculation) * 1.254 + *----------------------------------------------*/ 1.255 + 1.256 +typedef struct { 1.257 + pgl_point point; /* NOTE: point first to allow C cast to pgl_point */ 1.258 + int32 samples; 1.259 +} pgl_point_sc; 1.260 + 1.261 + 1.262 +/*----------------------------------------* 1.263 + * C functions on geographic data types * 1.264 + *----------------------------------------*/ 1.265 + 1.266 +/* round latitude or longitude to 12 digits after decimal point */ 1.267 +static inline double pgl_round(double val) { 1.268 + return round(val * 1e12) / 1e12; 1.269 +} 1.270 + 1.271 +/* compare two points */ 1.272 +/* (equality when same point on earth is described, otherwise an arbitrary 1.273 + linear order) */ 1.274 +static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) { 1.275 + double lon1, lon2; /* modified longitudes for special cases */ 1.276 + /* use latitude as first ordering criterion */ 1.277 + if (point1->lat < point2->lat) return -1; 1.278 + if (point1->lat > point2->lat) return 1; 1.279 + /* determine modified longitudes (considering special case of poles and 1.280 + 180th meridian which can be described as W180 or E180) */ 1.281 + if (point1->lat == -90 || point1->lat == 90) lon1 = 0; 1.282 + else if (point1->lon == 180) lon1 = -180; 1.283 + else lon1 = point1->lon; 1.284 + if (point2->lat == -90 || point2->lat == 90) lon2 = 0; 1.285 + else if (point2->lon == 180) lon2 = -180; 1.286 + else lon2 = point2->lon; 1.287 + /* use (modified) longitude as secondary ordering criterion */ 1.288 + if (lon1 < lon2) return -1; 1.289 + if (lon1 > lon2) return 1; 1.290 + /* no difference found, points are equal */ 1.291 + return 0; 1.292 +} 1.293 + 1.294 +/* compare two boxes */ 1.295 +/* (equality when same box on earth is described, otherwise an arbitrary linear 1.296 + order) */ 1.297 +static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) { 1.298 + /* two empty boxes are equal, and an empty box is always considered "less 1.299 + than" a non-empty box */ 1.300 + if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1; 1.301 + if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0; 1.302 + if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1; 1.303 + /* use southern border as first ordering criterion */ 1.304 + if (box1->lat_min < box2->lat_min) return -1; 1.305 + if (box1->lat_min > box2->lat_min) return 1; 1.306 + /* use northern border as second ordering criterion */ 1.307 + if (box1->lat_max < box2->lat_max) return -1; 1.308 + if (box1->lat_max > box2->lat_max) return 1; 1.309 + /* use western border as third ordering criterion */ 1.310 + if (box1->lon_min < box2->lon_min) return -1; 1.311 + if (box1->lon_min > box2->lon_min) return 1; 1.312 + /* use eastern border as fourth ordering criterion */ 1.313 + if (box1->lon_max < box2->lon_max) return -1; 1.314 + if (box1->lon_max > box2->lon_max) return 1; 1.315 + /* no difference found, boxes are equal */ 1.316 + return 0; 1.317 +} 1.318 + 1.319 +/* compare two circles */ 1.320 +/* (equality when same circle on earth is described, otherwise an arbitrary 1.321 + linear order) */ 1.322 +static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) { 1.323 + /* two circles with same infinite radius (positive or negative infinity) are 1.324 + considered equal independently of center point */ 1.325 + if ( 1.326 + !isfinite(circle1->radius) && !isfinite(circle2->radius) && 1.327 + circle1->radius == circle2->radius 1.328 + ) return 0; 1.329 + /* use radius as first ordering criterion */ 1.330 + if (circle1->radius < circle2->radius) return -1; 1.331 + if (circle1->radius > circle2->radius) return 1; 1.332 + /* use center point as secondary ordering criterion */ 1.333 + return pgl_point_cmp(&(circle1->center), &(circle2->center)); 1.334 +} 1.335 + 1.336 +/* set box to empty box*/ 1.337 +static void pgl_box_set_empty(pgl_box *box) { 1.338 + box->lat_min = INFINITY; 1.339 + box->lat_max = -INFINITY; 1.340 + box->lon_min = 0; 1.341 + box->lon_max = 0; 1.342 +} 1.343 + 1.344 +/* check if point is inside a box */ 1.345 +static bool pgl_point_in_box(pgl_point *point, pgl_box *box) { 1.346 + return ( 1.347 + point->lat >= box->lat_min && point->lat <= box->lat_max && ( 1.348 + (box->lon_min > box->lon_max) ? ( 1.349 + /* box crosses 180th meridian */ 1.350 + point->lon >= box->lon_min || point->lon <= box->lon_max 1.351 + ) : ( 1.352 + /* box does not cross the 180th meridian */ 1.353 + point->lon >= box->lon_min && point->lon <= box->lon_max 1.354 + ) 1.355 + ) 1.356 + ); 1.357 +} 1.358 + 1.359 +/* check if two boxes overlap */ 1.360 +static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) { 1.361 + return ( 1.362 + box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */ 1.363 + ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) && 1.364 + ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && ( 1.365 + ( 1.366 + /* check if one and only one box crosses the 180th meridian */ 1.367 + ((box1->lon_min > box1->lon_max) ? 1 : 0) ^ 1.368 + ((box2->lon_min > box2->lon_max) ? 1 : 0) 1.369 + ) ? ( 1.370 + /* exactly one box crosses the 180th meridian */ 1.371 + box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min || 1.372 + box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max 1.373 + ) : ( 1.374 + /* no box or both boxes cross the 180th meridian */ 1.375 + ( 1.376 + (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) && 1.377 + (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max) 1.378 + ) || 1.379 + /* handle W180 == E180 */ 1.380 + ( box1->lon_min == -180 && box2->lon_max == 180 ) || 1.381 + ( box2->lon_min == -180 && box1->lon_max == 180 ) 1.382 + ) 1.383 + ) 1.384 + ); 1.385 +} 1.386 + 1.387 +/* check unambiguousness of east/west orientation of cluster entries and set 1.388 + bounding circle of cluster */ 1.389 +static bool pgl_finalize_cluster(pgl_cluster *cluster) { 1.390 + int i, j; /* i: index of entry, j: index of point in entry */ 1.391 + int npoints; /* number of points in entry */ 1.392 + int total_npoints = 0; /* total number of points in cluster */ 1.393 + pgl_point *points; /* points in entry */ 1.394 + int lon_dir; /* first point of entry west (-1) or east (+1) */ 1.395 + double lon_break = 0; /* antipodal longitude of first point in entry */ 1.396 + double lon_min, lon_max; /* covered longitude range of entry */ 1.397 + double value; /* temporary variable */ 1.398 + /* reset bounding circle center to empty circle at 0/0 coordinates */ 1.399 + cluster->bounding.center.lat = 0; 1.400 + cluster->bounding.center.lon = 0; 1.401 + cluster->bounding.radius = -INFINITY; 1.402 + /* if cluster is not empty */ 1.403 + if (cluster->nentries != 0) { 1.404 + /* iterate over all cluster entries and ensure they each cover a longitude 1.405 + range less than 180 degrees */ 1.406 + for (i=0; i<cluster->nentries; i++) { 1.407 + /* get properties of entry */ 1.408 + npoints = cluster->entries[i].npoints; 1.409 + points = PGL_ENTRY_POINTS(cluster, i); 1.410 + /* get longitude of first point of entry */ 1.411 + value = points[0].lon; 1.412 + /* initialize lon_min and lon_max with longitude of first point */ 1.413 + lon_min = value; 1.414 + lon_max = value; 1.415 + /* determine east/west orientation of first point and calculate antipodal 1.416 + longitude (Note: rounding required here) */ 1.417 + if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); } 1.418 + else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); } 1.419 + else lon_dir = 0; 1.420 + /* iterate over all other points in entry */ 1.421 + for (j=1; j<npoints; j++) { 1.422 + /* consider longitude wrap-around */ 1.423 + value = points[j].lon; 1.424 + if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360); 1.425 + else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360); 1.426 + /* update lon_min and lon_max */ 1.427 + if (value < lon_min) lon_min = value; 1.428 + else if (value > lon_max) lon_max = value; 1.429 + /* return false if 180 degrees or more are covered */ 1.430 + if (lon_max - lon_min >= 180) return false; 1.431 + } 1.432 + } 1.433 + /* iterate over all points of all entries and calculate arbitrary center 1.434 + point for bounding circle (best if center point minimizes the radius, 1.435 + but some error is allowed here) */ 1.436 + for (i=0; i<cluster->nentries; i++) { 1.437 + /* get properties of entry */ 1.438 + npoints = cluster->entries[i].npoints; 1.439 + points = PGL_ENTRY_POINTS(cluster, i); 1.440 + /* check if first entry */ 1.441 + if (i==0) { 1.442 + /* get longitude of first point of first entry in whole cluster */ 1.443 + value = points[0].lon; 1.444 + /* initialize lon_min and lon_max with longitude of first point of 1.445 + first entry in whole cluster (used to determine if whole cluster 1.446 + covers a longitude range of 180 degrees or more) */ 1.447 + lon_min = value; 1.448 + lon_max = value; 1.449 + /* determine east/west orientation of first point and calculate 1.450 + antipodal longitude (Note: rounding not necessary here) */ 1.451 + if (value < 0) { lon_dir = -1; lon_break = value + 180; } 1.452 + else if (value > 0) { lon_dir = 1; lon_break = value - 180; } 1.453 + else lon_dir = 0; 1.454 + } 1.455 + /* iterate over all points in entry */ 1.456 + for (j=0; j<npoints; j++) { 1.457 + /* longitude wrap-around (Note: rounding not necessary here) */ 1.458 + value = points[j].lon; 1.459 + if (lon_dir < 0 && value > lon_break) value -= 360; 1.460 + else if (lon_dir > 0 && value < lon_break) value += 360; 1.461 + if (value < lon_min) lon_min = value; 1.462 + else if (value > lon_max) lon_max = value; 1.463 + /* set bounding circle to cover whole earth if 180 degrees or more are 1.464 + covered */ 1.465 + if (lon_max - lon_min >= 180) { 1.466 + cluster->bounding.center.lat = 0; 1.467 + cluster->bounding.center.lon = 0; 1.468 + cluster->bounding.radius = INFINITY; 1.469 + return true; 1.470 + } 1.471 + /* add point to bounding circle center (for average calculation) */ 1.472 + cluster->bounding.center.lat += points[j].lat; 1.473 + cluster->bounding.center.lon += value; 1.474 + } 1.475 + /* count total number of points */ 1.476 + total_npoints += npoints; 1.477 + } 1.478 + /* determine average latitude and longitude of cluster */ 1.479 + cluster->bounding.center.lat /= total_npoints; 1.480 + cluster->bounding.center.lon /= total_npoints; 1.481 + /* normalize longitude of center of cluster bounding circle */ 1.482 + if (cluster->bounding.center.lon < -180) { 1.483 + cluster->bounding.center.lon += 360; 1.484 + } 1.485 + else if (cluster->bounding.center.lon > 180) { 1.486 + cluster->bounding.center.lon -= 360; 1.487 + } 1.488 + /* round bounding circle center (useful if it is used by other functions) */ 1.489 + cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat); 1.490 + cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon); 1.491 + /* calculate radius of bounding circle */ 1.492 + for (i=0; i<cluster->nentries; i++) { 1.493 + npoints = cluster->entries[i].npoints; 1.494 + points = PGL_ENTRY_POINTS(cluster, i); 1.495 + for (j=0; j<npoints; j++) { 1.496 + value = pgl_distance( 1.497 + cluster->bounding.center.lat, cluster->bounding.center.lon, 1.498 + points[j].lat, points[j].lon 1.499 + ); 1.500 + if (value > cluster->bounding.radius) cluster->bounding.radius = value; 1.501 + } 1.502 + } 1.503 + } 1.504 + /* return true (east/west orientation is unambiguous) */ 1.505 + return true; 1.506 +} 1.507 + 1.508 +/* check if point is inside cluster */ 1.509 +/* (if point is on perimeter, then true is returned if and only if 1.510 + strict == false) */ 1.511 +static bool pgl_point_in_cluster( 1.512 + pgl_point *point, 1.513 + pgl_cluster *cluster, 1.514 + bool strict 1.515 +) { 1.516 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 1.517 + int entrytype; /* type of entry */ 1.518 + int npoints; /* number of points in entry */ 1.519 + pgl_point *points; /* array of points in entry */ 1.520 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 1.521 + double lon_break = 0; /* antipodal longitude of first vertex */ 1.522 + double lat0 = point->lat; /* latitude of point */ 1.523 + double lon0; /* (adjusted) longitude of point */ 1.524 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 1.525 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 1.526 + double lon; /* longitude of intersection */ 1.527 + int counter = 0; /* counter for intersections east of point */ 1.528 + /* iterate over all entries */ 1.529 + for (i=0; i<cluster->nentries; i++) { 1.530 + /* get type of entry */ 1.531 + entrytype = cluster->entries[i].entrytype; 1.532 + /* skip all entries but polygons if perimeters are excluded */ 1.533 + if (strict && entrytype != PGL_ENTRY_POLYGON) continue; 1.534 + /* get points of entry */ 1.535 + npoints = cluster->entries[i].npoints; 1.536 + points = PGL_ENTRY_POINTS(cluster, i); 1.537 + /* determine east/west orientation of first point of entry and calculate 1.538 + antipodal longitude */ 1.539 + lon_break = points[0].lon; 1.540 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 1.541 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 1.542 + else lon_dir = 0; 1.543 + /* get longitude of point */ 1.544 + lon0 = point->lon; 1.545 + /* consider longitude wrap-around for point */ 1.546 + if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360); 1.547 + else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360); 1.548 + /* iterate over all edges and vertices */ 1.549 + for (j=0; j<npoints; j++) { 1.550 + /* return if point is on vertex of polygon */ 1.551 + if (pgl_point_cmp(point, &(points[j])) == 0) return !strict; 1.552 + /* calculate index of next vertex */ 1.553 + k = (j+1) % npoints; 1.554 + /* skip last edge unless entry is (closed) outline or polygon */ 1.555 + if ( 1.556 + k == 0 && 1.557 + entrytype != PGL_ENTRY_OUTLINE && 1.558 + entrytype != PGL_ENTRY_POLYGON 1.559 + ) continue; 1.560 + /* use previously calculated values for lat1 and lon1 if possible */ 1.561 + if (j) { 1.562 + lat1 = lat2; 1.563 + lon1 = lon2; 1.564 + } else { 1.565 + /* otherwise get latitude and longitude values of first vertex */ 1.566 + lat1 = points[0].lat; 1.567 + lon1 = points[0].lon; 1.568 + /* and consider longitude wrap-around for first vertex */ 1.569 + if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360); 1.570 + else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360); 1.571 + } 1.572 + /* get latitude and longitude of next vertex */ 1.573 + lat2 = points[k].lat; 1.574 + lon2 = points[k].lon; 1.575 + /* consider longitude wrap-around for next vertex */ 1.576 + if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360); 1.577 + else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360); 1.578 + /* return if point is on horizontal (west to east) edge of polygon */ 1.579 + if ( 1.580 + lat0 == lat1 && lat0 == lat2 && 1.581 + ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) ) 1.582 + ) return !strict; 1.583 + /* check if edge crosses east/west line of point */ 1.584 + if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) { 1.585 + /* calculate longitude of intersection */ 1.586 + lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1); 1.587 + /* return if intersection goes (approximately) through point */ 1.588 + if (pgl_round(lon) == lon0) return !strict; 1.589 + /* count intersection if east of point and entry is polygon*/ 1.590 + if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++; 1.591 + } 1.592 + } 1.593 + } 1.594 + /* return true if number of intersections is odd */ 1.595 + return counter & 1; 1.596 +} 1.597 + 1.598 +/* check if all points of the second cluster are strictly inside the first 1.599 + cluster */ 1.600 +static inline bool pgl_all_cluster_points_strictly_in_cluster( 1.601 + pgl_cluster *outer, pgl_cluster *inner 1.602 +) { 1.603 + int i, j; /* i: entry, j: point in entry */ 1.604 + int npoints; /* number of points in entry */ 1.605 + pgl_point *points; /* array of points in entry */ 1.606 + /* iterate over all entries of "inner" cluster */ 1.607 + for (i=0; i<inner->nentries; i++) { 1.608 + /* get properties of entry */ 1.609 + npoints = inner->entries[i].npoints; 1.610 + points = PGL_ENTRY_POINTS(inner, i); 1.611 + /* iterate over all points in entry of "inner" cluster */ 1.612 + for (j=0; j<npoints; j++) { 1.613 + /* return false if one point of inner cluster is not in outer cluster */ 1.614 + if (!pgl_point_in_cluster(points+j, outer, true)) return false; 1.615 + } 1.616 + } 1.617 + /* otherwise return true */ 1.618 + return true; 1.619 +} 1.620 + 1.621 +/* check if any point the second cluster is inside the first cluster */ 1.622 +static inline bool pgl_any_cluster_points_in_cluster( 1.623 + pgl_cluster *outer, pgl_cluster *inner 1.624 +) { 1.625 + int i, j; /* i: entry, j: point in entry */ 1.626 + int npoints; /* number of points in entry */ 1.627 + pgl_point *points; /* array of points in entry */ 1.628 + /* iterate over all entries of "inner" cluster */ 1.629 + for (i=0; i<inner->nentries; i++) { 1.630 + /* get properties of entry */ 1.631 + npoints = inner->entries[i].npoints; 1.632 + points = PGL_ENTRY_POINTS(inner, i); 1.633 + /* iterate over all points in entry of "inner" cluster */ 1.634 + for (j=0; j<npoints; j++) { 1.635 + /* return true if one point of inner cluster is in outer cluster */ 1.636 + if (pgl_point_in_cluster(points+j, outer, false)) return true; 1.637 + } 1.638 + } 1.639 + /* otherwise return false */ 1.640 + return false; 1.641 +} 1.642 + 1.643 +/* check if line segment strictly crosses line (not just touching) */ 1.644 +static inline bool pgl_lseg_crosses_line( 1.645 + double seg_x1, double seg_y1, double seg_x2, double seg_y2, 1.646 + double line_x1, double line_y1, double line_x2, double line_y2 1.647 +) { 1.648 + return ( 1.649 + ( 1.650 + (seg_x1-line_x1) * (line_y2-line_y1) - 1.651 + (seg_y1-line_y1) * (line_x2-line_x1) 1.652 + ) * ( 1.653 + (seg_x2-line_x1) * (line_y2-line_y1) - 1.654 + (seg_y2-line_y1) * (line_x2-line_x1) 1.655 + ) 1.656 + ) < 0; 1.657 +} 1.658 + 1.659 +/* check if paths and outlines of two clusters strictly overlap (not just 1.660 + touching) */ 1.661 +static bool pgl_outlines_overlap( 1.662 + pgl_cluster *cluster1, pgl_cluster *cluster2 1.663 +) { 1.664 + int i1, j1, k1; /* i: entry, j: point in entry, k: next point in entry */ 1.665 + int i2, j2, k2; 1.666 + int entrytype1, entrytype2; /* type of entry */ 1.667 + int npoints1, npoints2; /* number of points in entry */ 1.668 + pgl_point *points1; /* array of points in entry of cluster1 */ 1.669 + pgl_point *points2; /* array of points in entry of cluster2 */ 1.670 + int lon_dir1, lon_dir2; /* first vertex west (-1) or east (+1) */ 1.671 + double lon_break1, lon_break2; /* antipodal longitude of first vertex */ 1.672 + double lat11, lon11; /* latitude and (adjusted) longitude of vertex */ 1.673 + double lat12, lon12; /* latitude and (adjusted) longitude of next vertex */ 1.674 + double lat21, lon21; /* latitude and (adjusted) longitudes for cluster2 */ 1.675 + double lat22, lon22; 1.676 + double wrapvalue; /* temporary helper value to adjust wrap-around */ 1.677 + /* iterate over all entries of cluster1 */ 1.678 + for (i1=0; i1<cluster1->nentries; i1++) { 1.679 + /* get properties of entry in cluster1 and skip points */ 1.680 + npoints1 = cluster1->entries[i1].npoints; 1.681 + if (npoints1 < 2) continue; 1.682 + entrytype1 = cluster1->entries[i1].entrytype; 1.683 + points1 = PGL_ENTRY_POINTS(cluster1, i1); 1.684 + /* determine east/west orientation of first point and calculate antipodal 1.685 + longitude */ 1.686 + lon_break1 = points1[0].lon; 1.687 + if (lon_break1 < 0) { 1.688 + lon_dir1 = -1; 1.689 + lon_break1 = pgl_round(lon_break1 + 180); 1.690 + } else if (lon_break1 > 0) { 1.691 + lon_dir1 = 1; 1.692 + lon_break1 = pgl_round(lon_break1 - 180); 1.693 + } else lon_dir1 = 0; 1.694 + /* iterate over all edges and vertices in cluster1 */ 1.695 + for (j1=0; j1<npoints1; j1++) { 1.696 + /* calculate index of next vertex */ 1.697 + k1 = (j1+1) % npoints1; 1.698 + /* skip last edge unless entry is (closed) outline or polygon */ 1.699 + if ( 1.700 + k1 == 0 && 1.701 + entrytype1 != PGL_ENTRY_OUTLINE && 1.702 + entrytype1 != PGL_ENTRY_POLYGON 1.703 + ) continue; 1.704 + /* use previously calculated values for lat1 and lon1 if possible */ 1.705 + if (j1) { 1.706 + lat11 = lat12; 1.707 + lon11 = lon12; 1.708 + } else { 1.709 + /* otherwise get latitude and longitude values of first vertex */ 1.710 + lat11 = points1[0].lat; 1.711 + lon11 = points1[0].lon; 1.712 + /* and consider longitude wrap-around for first vertex */ 1.713 + if (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360); 1.714 + else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360); 1.715 + } 1.716 + /* get latitude and longitude of next vertex */ 1.717 + lat12 = points1[k1].lat; 1.718 + lon12 = points1[k1].lon; 1.719 + /* consider longitude wrap-around for next vertex */ 1.720 + if (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360); 1.721 + else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360); 1.722 + /* skip degenerated edges */ 1.723 + if (lat11 == lat12 && lon11 == lon12) continue; 1.724 + /* iterate over all entries of cluster2 */ 1.725 + for (i2=0; i2<cluster2->nentries; i2++) { 1.726 + /* get points and number of points of entry in cluster2 */ 1.727 + npoints2 = cluster2->entries[i2].npoints; 1.728 + if (npoints2 < 2) continue; 1.729 + entrytype2 = cluster2->entries[i2].entrytype; 1.730 + points2 = PGL_ENTRY_POINTS(cluster2, i2); 1.731 + /* determine east/west orientation of first point and calculate antipodal 1.732 + longitude */ 1.733 + lon_break2 = points2[0].lon; 1.734 + if (lon_break2 < 0) { 1.735 + lon_dir2 = -1; 1.736 + lon_break2 = pgl_round(lon_break2 + 180); 1.737 + } else if (lon_break2 > 0) { 1.738 + lon_dir2 = 1; 1.739 + lon_break2 = pgl_round(lon_break2 - 180); 1.740 + } else lon_dir2 = 0; 1.741 + /* iterate over all edges and vertices in cluster2 */ 1.742 + for (j2=0; j2<npoints2; j2++) { 1.743 + /* calculate index of next vertex */ 1.744 + k2 = (j2+1) % npoints2; 1.745 + /* skip last edge unless entry is (closed) outline or polygon */ 1.746 + if ( 1.747 + k2 == 0 && 1.748 + entrytype2 != PGL_ENTRY_OUTLINE && 1.749 + entrytype2 != PGL_ENTRY_POLYGON 1.750 + ) continue; 1.751 + /* use previously calculated values for lat1 and lon1 if possible */ 1.752 + if (j2) { 1.753 + lat21 = lat22; 1.754 + lon21 = lon22; 1.755 + } else { 1.756 + /* otherwise get latitude and longitude values of first vertex */ 1.757 + lat21 = points2[0].lat; 1.758 + lon21 = points2[0].lon; 1.759 + /* and consider longitude wrap-around for first vertex */ 1.760 + if (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360); 1.761 + else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360); 1.762 + } 1.763 + /* get latitude and longitude of next vertex */ 1.764 + lat22 = points2[k2].lat; 1.765 + lon22 = points2[k2].lon; 1.766 + /* consider longitude wrap-around for next vertex */ 1.767 + if (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360); 1.768 + else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360); 1.769 + /* skip degenerated edges */ 1.770 + if (lat21 == lat22 && lon21 == lon22) continue; 1.771 + /* perform another wrap-around where necessary */ 1.772 + /* TODO: improve performance of whole wrap-around mechanism */ 1.773 + wrapvalue = (lon21 + lon22) - (lon11 + lon12); 1.774 + if (wrapvalue > 360) { 1.775 + lon21 = pgl_round(lon21 - 360); 1.776 + lon22 = pgl_round(lon22 - 360); 1.777 + } else if (wrapvalue < -360) { 1.778 + lon21 = pgl_round(lon21 + 360); 1.779 + lon22 = pgl_round(lon22 + 360); 1.780 + } 1.781 + /* return true if segments overlap */ 1.782 + if ( 1.783 + pgl_lseg_crosses_line( 1.784 + lat11, lon11, lat12, lon12, 1.785 + lat21, lon21, lat22, lon22 1.786 + ) && pgl_lseg_crosses_line( 1.787 + lat21, lon21, lat22, lon22, 1.788 + lat11, lon11, lat12, lon12 1.789 + ) 1.790 + ) { 1.791 + return true; 1.792 + } 1.793 + } 1.794 + } 1.795 + } 1.796 + } 1.797 + /* otherwise return false */ 1.798 + return false; 1.799 +} 1.800 + 1.801 +/* check if second cluster is completely contained in first cluster */ 1.802 +static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) { 1.803 + if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false; 1.804 + if (pgl_any_cluster_points_in_cluster(inner, outer)) return false; 1.805 + if (pgl_outlines_overlap(outer, inner)) return false; 1.806 + return true; 1.807 +} 1.808 + 1.809 +/* check if two clusters overlap */ 1.810 +static bool pgl_clusters_overlap( 1.811 + pgl_cluster *cluster1, pgl_cluster *cluster2 1.812 +) { 1.813 + if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true; 1.814 + if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true; 1.815 + if (pgl_outlines_overlap(cluster1, cluster2)) return true; 1.816 + return false; 1.817 +} 1.818 + 1.819 +/* calculate (approximate) distance between point and cluster */ 1.820 +static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) { 1.821 + double comp; /* square of compression of meridians */ 1.822 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 1.823 + int entrytype; /* type of entry */ 1.824 + int npoints; /* number of points in entry */ 1.825 + pgl_point *points; /* array of points in entry */ 1.826 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 1.827 + double lon_break = 0; /* antipodal longitude of first vertex */ 1.828 + double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */ 1.829 + double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */ 1.830 + double lat0 = point->lat; /* latitude of point */ 1.831 + double lon0; /* (adjusted) longitude of point */ 1.832 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 1.833 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 1.834 + double s; /* scalar for vector calculations */ 1.835 + double dist; /* distance calculated in one step */ 1.836 + double min_dist = INFINITY; /* minimum distance */ 1.837 + /* distance is zero if point is contained in cluster */ 1.838 + if (pgl_point_in_cluster(point, cluster, false)) return 0; 1.839 + /* calculate approximate square compression of meridians */ 1.840 + comp = cos((lat0 / 180.0) * M_PI); 1.841 + comp *= comp; 1.842 + /* calculate exact square compression of meridians */ 1.843 + comp *= ( 1.844 + (1.0 - PGL_EPS2 * (1.0-comp)) * 1.845 + (1.0 - PGL_EPS2 * (1.0-comp)) / 1.846 + (PGL_SUBEPS2 * PGL_SUBEPS2) 1.847 + ); 1.848 + /* iterate over all entries */ 1.849 + for (i=0; i<cluster->nentries; i++) { 1.850 + /* get properties of entry */ 1.851 + entrytype = cluster->entries[i].entrytype; 1.852 + npoints = cluster->entries[i].npoints; 1.853 + points = PGL_ENTRY_POINTS(cluster, i); 1.854 + /* determine east/west orientation of first point of entry and calculate 1.855 + antipodal longitude */ 1.856 + lon_break = points[0].lon; 1.857 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 1.858 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 1.859 + else lon_dir = 0; 1.860 + /* determine covered longitude range */ 1.861 + for (j=0; j<npoints; j++) { 1.862 + /* get longitude of vertex */ 1.863 + lon1 = points[j].lon; 1.864 + /* adjust longitude to fix potential wrap-around */ 1.865 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 1.866 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 1.867 + /* update minimum and maximum longitude of polygon */ 1.868 + if (j == 0 || lon1 < lon_min) lon_min = lon1; 1.869 + if (j == 0 || lon1 > lon_max) lon_max = lon1; 1.870 + } 1.871 + /* adjust longitude wrap-around according to full longitude range */ 1.872 + lon_break = (lon_max + lon_min) / 2; 1.873 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 1.874 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 1.875 + /* get longitude of point */ 1.876 + lon0 = point->lon; 1.877 + /* consider longitude wrap-around for point */ 1.878 + if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360; 1.879 + else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360; 1.880 + /* iterate over all edges and vertices */ 1.881 + for (j=0; j<npoints; j++) { 1.882 + /* use previously calculated values for lat1 and lon1 if possible */ 1.883 + if (j) { 1.884 + lat1 = lat2; 1.885 + lon1 = lon2; 1.886 + } else { 1.887 + /* otherwise get latitude and longitude values of first vertex */ 1.888 + lat1 = points[0].lat; 1.889 + lon1 = points[0].lon; 1.890 + /* and consider longitude wrap-around for first vertex */ 1.891 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 1.892 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 1.893 + } 1.894 + /* calculate distance to vertex */ 1.895 + dist = pgl_distance(lat0, lon0, lat1, lon1); 1.896 + /* store calculated distance if smallest */ 1.897 + if (dist < min_dist) min_dist = dist; 1.898 + /* calculate index of next vertex */ 1.899 + k = (j+1) % npoints; 1.900 + /* skip last edge unless entry is (closed) outline or polygon */ 1.901 + if ( 1.902 + k == 0 && 1.903 + entrytype != PGL_ENTRY_OUTLINE && 1.904 + entrytype != PGL_ENTRY_POLYGON 1.905 + ) continue; 1.906 + /* get latitude and longitude of next vertex */ 1.907 + lat2 = points[k].lat; 1.908 + lon2 = points[k].lon; 1.909 + /* consider longitude wrap-around for next vertex */ 1.910 + if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360; 1.911 + else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360; 1.912 + /* go to next vertex and edge if edge is degenerated */ 1.913 + if (lat1 == lat2 && lon1 == lon2) continue; 1.914 + /* otherwise test if point can be projected onto edge of polygon */ 1.915 + s = ( 1.916 + ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) / 1.917 + ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1)) 1.918 + ); 1.919 + /* go to next vertex and edge if point cannot be projected */ 1.920 + if (!(s > 0 && s < 1)) continue; 1.921 + /* calculate distance from original point to projected point */ 1.922 + dist = pgl_distance( 1.923 + lat0, lon0, 1.924 + lat1 + s * (lat2-lat1), 1.925 + lon1 + s * (lon2-lon1) 1.926 + ); 1.927 + /* store calculated distance if smallest */ 1.928 + if (dist < min_dist) min_dist = dist; 1.929 + } 1.930 + } 1.931 + /* return minimum distance */ 1.932 + return min_dist; 1.933 +} 1.934 + 1.935 +/* calculate (approximate) distance between two clusters */ 1.936 +static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) { 1.937 + int i, j; /* i: entry, j: point in entry */ 1.938 + int npoints; /* number of points in entry */ 1.939 + pgl_point *points; /* array of points in entry */ 1.940 + double dist; /* distance calculated in one step */ 1.941 + double min_dist = INFINITY; /* minimum distance */ 1.942 + /* consider distance from each point in one cluster to the whole other */ 1.943 + for (i=0; i<cluster1->nentries; i++) { 1.944 + npoints = cluster1->entries[i].npoints; 1.945 + points = PGL_ENTRY_POINTS(cluster1, i); 1.946 + for (j=0; j<npoints; j++) { 1.947 + dist = pgl_point_cluster_distance(points+j, cluster2); 1.948 + if (dist == 0) return dist; 1.949 + if (dist < min_dist) min_dist = dist; 1.950 + } 1.951 + } 1.952 + /* consider distance from each point in other cluster to the first cluster */ 1.953 + for (i=0; i<cluster2->nentries; i++) { 1.954 + npoints = cluster2->entries[i].npoints; 1.955 + points = PGL_ENTRY_POINTS(cluster2, i); 1.956 + for (j=0; j<npoints; j++) { 1.957 + dist = pgl_point_cluster_distance(points+j, cluster1); 1.958 + if (dist == 0) return dist; 1.959 + if (dist < min_dist) min_dist = dist; 1.960 + } 1.961 + } 1.962 + return min_dist; 1.963 +} 1.964 + 1.965 +/* estimator function for distance between box and point */ 1.966 +/* always returns a smaller value than actually correct or zero */ 1.967 +static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) { 1.968 + double dlon; /* longitude range of box (delta longitude) */ 1.969 + double distance; /* return value */ 1.970 + /* return infinity if box is empty */ 1.971 + if (box->lat_min > box->lat_max) return INFINITY; 1.972 + /* return zero if point is inside box */ 1.973 + if (pgl_point_in_box(point, box)) return 0; 1.974 + /* calculate delta longitude */ 1.975 + dlon = box->lon_max - box->lon_min; 1.976 + if (dlon < 0) dlon += 360; /* 180th meridian crossed */ 1.977 + /* if delta longitude is greater than 150 degrees, perform safe fall-back */ 1.978 + if (dlon > 150) return 0; 1.979 + /* calculate lower limit for distance (formula below requires dlon <= 150) */ 1.980 + /* TODO: provide better estimation function to improve performance */ 1.981 + distance = ( 1.982 + (1.0-PGL_SPHEROID_F) * /* safety margin due to flattening and approx. */ 1.983 + pgl_distance( 1.984 + point->lat, 1.985 + point->lon, 1.986 + (box->lat_min + box->lat_max) / 2, 1.987 + box->lon_min + dlon/2 1.988 + ) 1.989 + ) - pgl_distance( 1.990 + box->lat_min, box->lon_min, 1.991 + box->lat_max, box->lon_max 1.992 + ); 1.993 + /* truncate negative results to zero */ 1.994 + if (distance <= 0) distance = 0; 1.995 + /* return result */ 1.996 + return distance; 1.997 +} 1.998 + 1.999 + 1.1000 +/*------------------------------------------------------------* 1.1001 + * Functions using numerical integration (Monte Carlo like) * 1.1002 + *------------------------------------------------------------*/ 1.1003 + 1.1004 +/* half of (spherical) earth's surface area */ 1.1005 +#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI) 1.1006 + 1.1007 +/* golden angle in radians */ 1.1008 +#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0)) 1.1009 + 1.1010 +/* create a list of sample points covering a bounding circle 1.1011 + and return covered area */ 1.1012 +static double pgl_sample_points( 1.1013 + pgl_point *center, /* center of bounding circle */ 1.1014 + double radius, /* radius of bounding circle */ 1.1015 + int samples, /* number of sample points (MUST be positive!) */ 1.1016 + pgl_point *result /* pointer to result array */ 1.1017 +) { 1.1018 + double double_share = 2.0; /* double of covered share of earth's surface */ 1.1019 + double double_share_div_samples; /* double_share divided by sample count */ 1.1020 + int i; 1.1021 + double t; /* parameter of spiral laid on (spherical) earth's surface */ 1.1022 + double x, y, z; /* normalized coordinates of point on non-rotated spiral */ 1.1023 + double sin_phi; /* sine of sph. coordinate of point of non-rotated spiral */ 1.1024 + double lambda; /* other sph. coordinate of point of non-rotated spiral */ 1.1025 + double rot = (0.5 - center->lat / 180.0) * M_PI; /* needed rot. (in rad) */ 1.1026 + double cos_rot = cos(rot); /* cosine of rotation by latitude */ 1.1027 + double sin_rot = sin(rot); /* sine of rotation by latitude */ 1.1028 + double x_rot, z_rot; /* normalized coordinates of point on rotated spiral */ 1.1029 + double center_lon = center->lon; /* second rotation in degree */ 1.1030 + /* add safety margin to bounding circle because of spherical approximation */ 1.1031 + radius *= PGL_SPHEROID_A / PGL_RADIUS; 1.1032 + /* if whole earth is covered, use initialized value, otherwise calculate 1.1033 + share of covered area (multiplied by 2) */ 1.1034 + if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS); 1.1035 + /* divide double_share by sample count for later calculations */ 1.1036 + double_share_div_samples = double_share / samples; 1.1037 + /* generate sample points */ 1.1038 + for (i=0; i<samples; i++) { 1.1039 + /* use an offset of 1/2 to avoid too dense clustering at spiral center */ 1.1040 + t = 0.5 + i; 1.1041 + /* calculate normalized coordinates of point on non-rotated spiral */ 1.1042 + z = 1.0 - double_share_div_samples * t; 1.1043 + sin_phi = sqrt(1.0 - z*z); 1.1044 + lambda = t * PGL_GOLDEN_ANGLE; 1.1045 + x = sin_phi * cos(lambda); 1.1046 + y = sin_phi * sin(lambda); 1.1047 + /* rotate spiral by latitude value of bounding circle */ 1.1048 + x_rot = cos_rot * x + sin_rot * z; 1.1049 + z_rot = cos_rot * z - sin_rot * x; 1.1050 + /* set resulting sample point in result array */ 1.1051 + /* (while performing second rotation by bounding circle longitude) */ 1.1052 + result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI); 1.1053 + result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI); 1.1054 + } 1.1055 + /* return covered area */ 1.1056 + return PGL_HALF_SURFACE * double_share; 1.1057 +} 1.1058 + 1.1059 +/* fair distance between point and cluster (see README file for explanation) */ 1.1060 +/* NOTE: sample count passed as third argument MUST be positive */ 1.1061 +static double pgl_fair_distance( 1.1062 + pgl_point *point, pgl_cluster *cluster, int samples 1.1063 +) { 1.1064 + double distance; /* shortest distance from point to cluster */ 1.1065 + pgl_point *points; /* sample points for numerical integration */ 1.1066 + double area; /* area covered by sample points */ 1.1067 + int i; 1.1068 + int inner = 0; /* number of sample points within cluster */ 1.1069 + int outer = 0; /* number of sample points outside cluster but 1.1070 + within cluster enlarged by distance */ 1.1071 + double result; 1.1072 + /* calculate shortest distance from point to cluster */ 1.1073 + distance = pgl_point_cluster_distance(point, cluster); 1.1074 + /* if cluster consists of a single point or has no bounding circle with 1.1075 + positive radius, simply return distance */ 1.1076 + if ( 1.1077 + (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) || 1.1078 + !(cluster->bounding.radius > 0) 1.1079 + ) return distance; 1.1080 + /* if cluster consists of two points which are twice as far apart, return 1.1081 + distance between point and cluster multiplied by square root of two */ 1.1082 + if ( 1.1083 + cluster->nentries == 2 && 1.1084 + cluster->entries[0].entrytype == PGL_ENTRY_POINT && 1.1085 + cluster->entries[1].entrytype == PGL_ENTRY_POINT && 1.1086 + pgl_distance( 1.1087 + PGL_ENTRY_POINTS(cluster, 0)[0].lat, 1.1088 + PGL_ENTRY_POINTS(cluster, 0)[0].lon, 1.1089 + PGL_ENTRY_POINTS(cluster, 1)[0].lat, 1.1090 + PGL_ENTRY_POINTS(cluster, 1)[0].lon 1.1091 + ) >= 2.0 * distance 1.1092 + ) { 1.1093 + return distance * M_SQRT2; 1.1094 + } 1.1095 + /* otherwise create sample points for numerical integration and determine 1.1096 + area covered by sample points */ 1.1097 + points = palloc(samples * sizeof(pgl_point)); 1.1098 + area = pgl_sample_points( 1.1099 + &cluster->bounding.center, 1.1100 + cluster->bounding.radius + distance, /* pad bounding circle by distance */ 1.1101 + samples, 1.1102 + points 1.1103 + ); 1.1104 + /* perform numerical integration */ 1.1105 + if (distance > 0) { 1.1106 + /* point (that was passed as argument) is outside cluster */ 1.1107 + for (i=0; i<samples; i++) { 1.1108 + /* count sample points within cluster */ 1.1109 + if (pgl_point_in_cluster(points+i, cluster, true)) inner++; 1.1110 + /* count sample points outside of cluster but within cluster enlarged by 1.1111 + distance between point (that was passed as argument) and cluster */ 1.1112 + else if ( 1.1113 + pgl_point_cluster_distance(points+i, cluster) < distance 1.1114 + ) outer++; 1.1115 + } 1.1116 + } else { 1.1117 + /* if point is within cluster, just count sample points within cluster */ 1.1118 + for (i=0; i<samples; i++) { 1.1119 + if (pgl_point_in_cluster(points+i, cluster, true)) inner++; 1.1120 + } 1.1121 + } 1.1122 + /* release memory for sample points needed for numerical integration */ 1.1123 + pfree(points); 1.1124 + /* if enlargement was less than doubling the area, then combine inner and 1.1125 + outer sample point counts with different weighting */ 1.1126 + /* (ensures fairness in such a way that the integral of the squared result 1.1127 + over all possible point parameters is independent of the cluster) */ 1.1128 + if (outer < inner) result = (2*inner + 4*outer) / 3.0; 1.1129 + /* otherwise weigh inner and outer points the same */ 1.1130 + else result = inner + outer; 1.1131 + /* convert area into distance (i.e. radius of a circle with the same area) */ 1.1132 + result = sqrt(area * (result / samples) / M_PI); 1.1133 + /* return result only if it is greater than the distance between point and 1.1134 + cluster to avoid unexpected results because of errors due to limited 1.1135 + precision */ 1.1136 + if (result > distance) return result; 1.1137 + /* otherwise return distance between point and cluster */ 1.1138 + else return distance; 1.1139 +} 1.1140 + 1.1141 + 1.1142 +/*-------------------------------------------------* 1.1143 + * geographic index based on space-filling curve * 1.1144 + *-------------------------------------------------*/ 1.1145 + 1.1146 +/* number of bytes used for geographic (center) position in keys */ 1.1147 +#define PGL_KEY_LATLON_BYTELEN 7 1.1148 + 1.1149 +/* maximum reference value for logarithmic size of geographic objects */ 1.1150 +#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */ 1.1151 + 1.1152 +/* pointer to index key (either pgl_pointkey or pgl_areakey) */ 1.1153 +typedef unsigned char *pgl_keyptr; 1.1154 + 1.1155 +/* index key for points (objects with zero area) on the spheroid */ 1.1156 +/* bit 0..55: interspersed bits of latitude and longitude, 1.1157 + bit 56..57: always zero, 1.1158 + bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */ 1.1159 +typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1]; 1.1160 + 1.1161 +/* index key for geographic objects on spheroid with area greater than zero */ 1.1162 +/* bit 0..55: interspersed bits of latitude and longitude of center point, 1.1163 + bit 56: always set to 1, 1.1164 + bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.), 1.1165 + bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to 1.1166 + PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth 1.1167 + = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL 1.1168 + (with interspersed bits = 0 and node depth = 0) for keys which 1.1169 + cover both empty and non-empty objects */ 1.1170 + 1.1171 +typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2]; 1.1172 + 1.1173 +/* helper macros for reading/writing index keys */ 1.1174 +#define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN 1.1175 +#define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1) 1.1176 +#define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8) 1.1177 +#define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1) 1.1178 +#define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1) 1.1179 +#define PGL_AREAKEY_TYPEMASK 0x80 1.1180 +#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8))) 1.1181 +#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \ 1.1182 + ( PGL_KEY_LATLONBIT(key1, n) ^ \ 1.1183 + PGL_KEY_LATLONBIT(key2, n) ) 1.1184 +#define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 1.1185 + PGL_AREAKEY_TYPEMASK) 1.1186 +#define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 1.1187 + (PGL_AREAKEY_TYPEMASK-1)) 1.1188 +#define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET]) 1.1189 +#define PGL_KEY_OBJSIZE_EMPTY 126 1.1190 +#define PGL_KEY_OBJSIZE_UNIVERSAL 127 1.1191 +#define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \ 1.1192 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 1.1193 + PGL_KEY_OBJSIZE_EMPTY ) 1.1194 +#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \ 1.1195 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 1.1196 + PGL_KEY_OBJSIZE_UNIVERSAL ) 1.1197 + 1.1198 +/* set area key to match empty objects only */ 1.1199 +static void pgl_key_set_empty(pgl_keyptr key) { 1.1200 + memset(key, 0, sizeof(pgl_areakey)); 1.1201 + /* Note: setting node depth to maximum is required for picksplit function */ 1.1202 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 1.1203 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY; 1.1204 +} 1.1205 + 1.1206 +/* set area key to match any object (including empty objects) */ 1.1207 +static void pgl_key_set_universal(pgl_keyptr key) { 1.1208 + memset(key, 0, sizeof(pgl_areakey)); 1.1209 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK; 1.1210 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL; 1.1211 +} 1.1212 + 1.1213 +/* convert a point on earth into a max-depth key to be used in index */ 1.1214 +static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) { 1.1215 + double lat = point->lat; 1.1216 + double lon = point->lon; 1.1217 + int i; 1.1218 + /* clear latitude and longitude bits */ 1.1219 + memset(key, 0, PGL_KEY_LATLON_BYTELEN); 1.1220 + /* set node depth to maximum and type bit to zero */ 1.1221 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH; 1.1222 + /* iterate over all latitude/longitude bit pairs */ 1.1223 + for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) { 1.1224 + /* determine latitude bit */ 1.1225 + if (lat >= 0) { 1.1226 + key[i/4] |= 0x80 >> (2*(i%4)); 1.1227 + lat *= 2; lat -= 90; 1.1228 + } else { 1.1229 + lat *= 2; lat += 90; 1.1230 + } 1.1231 + /* determine longitude bit */ 1.1232 + if (lon >= 0) { 1.1233 + key[i/4] |= 0x80 >> (2*(i%4)+1); 1.1234 + lon *= 2; lon -= 180; 1.1235 + } else { 1.1236 + lon *= 2; lon += 180; 1.1237 + } 1.1238 + } 1.1239 +} 1.1240 + 1.1241 +/* convert a circle on earth into a max-depth key to be used in an index */ 1.1242 +static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) { 1.1243 + /* handle special case of empty circle */ 1.1244 + if (circle->radius < 0) { 1.1245 + pgl_key_set_empty(key); 1.1246 + return; 1.1247 + } 1.1248 + /* perform same action as for point keys */ 1.1249 + pgl_point_to_key(&(circle->center), key); 1.1250 + /* but overwrite type and node depth to fit area index key */ 1.1251 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 1.1252 + /* check if radius is greater than (or equal to) reference size */ 1.1253 + /* (treat equal values as greater values for numerical safety) */ 1.1254 + if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) { 1.1255 + /* if yes, set logarithmic size to zero */ 1.1256 + key[PGL_KEY_OBJSIZE_OFFSET] = 0; 1.1257 + } else { 1.1258 + /* otherwise, determine logarithmic size iteratively */ 1.1259 + /* (one step is equivalent to a factor of sqrt(2)) */ 1.1260 + double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2; 1.1261 + int objsize = 1; 1.1262 + while (objsize < PGL_AREAKEY_MAXOBJSIZE) { 1.1263 + /* stop when radius is greater than (or equal to) adjusted reference */ 1.1264 + /* (treat equal values as greater values for numerical safety) */ 1.1265 + if (circle->radius >= reference) break; 1.1266 + reference /= M_SQRT2; 1.1267 + objsize++; 1.1268 + } 1.1269 + /* set logarithmic size to determined value */ 1.1270 + key[PGL_KEY_OBJSIZE_OFFSET] = objsize; 1.1271 + } 1.1272 +} 1.1273 + 1.1274 +/* check if one key is subkey of another key or vice versa */ 1.1275 +static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) { 1.1276 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 1.1277 + /* determine smallest depth */ 1.1278 + int depth1 = PGL_KEY_NODEDEPTH(key1); 1.1279 + int depth2 = PGL_KEY_NODEDEPTH(key2); 1.1280 + int depth = (depth1 < depth2) ? depth1 : depth2; 1.1281 + /* check if keys are area keys (assuming that both keys have same type) */ 1.1282 + if (PGL_KEY_IS_AREAKEY(key1)) { 1.1283 + int j = 0; /* bit offset for logarithmic object size bits */ 1.1284 + int k = 0; /* bit offset for latitude and longitude */ 1.1285 + /* fetch logarithmic object size information */ 1.1286 + int objsize1 = PGL_KEY_OBJSIZE(key1); 1.1287 + int objsize2 = PGL_KEY_OBJSIZE(key2); 1.1288 + /* handle special cases for empty objects (universal and empty keys) */ 1.1289 + if ( 1.1290 + objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL || 1.1291 + objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL 1.1292 + ) return true; 1.1293 + if ( 1.1294 + objsize1 == PGL_KEY_OBJSIZE_EMPTY || 1.1295 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 1.1296 + ) return objsize1 == objsize2; 1.1297 + /* iterate through key bits */ 1.1298 + for (i=0; i<depth; i++) { 1.1299 + /* every second bit is a bit describing the object size */ 1.1300 + if (i%2 == 0) { 1.1301 + /* check if object size bit is different in both keys (objsize1 and 1.1302 + objsize2 describe the minimum index when object size bit is set) */ 1.1303 + if ( 1.1304 + (objsize1 <= j && objsize2 > j) || 1.1305 + (objsize2 <= j && objsize1 > j) 1.1306 + ) { 1.1307 + /* bit differs, therefore keys are in separate branches */ 1.1308 + return false; 1.1309 + } 1.1310 + /* increase bit counter for object size bits */ 1.1311 + j++; 1.1312 + } 1.1313 + /* all other bits describe latitude and longitude */ 1.1314 + else { 1.1315 + /* check if bit differs in both keys */ 1.1316 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) { 1.1317 + /* bit differs, therefore keys are in separate branches */ 1.1318 + return false; 1.1319 + } 1.1320 + /* increase bit counter for latitude/longitude bits */ 1.1321 + k++; 1.1322 + } 1.1323 + } 1.1324 + } 1.1325 + /* if not, keys are point keys */ 1.1326 + else { 1.1327 + /* iterate through key bits */ 1.1328 + for (i=0; i<depth; i++) { 1.1329 + /* check if bit differs in both keys */ 1.1330 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) { 1.1331 + /* bit differs, therefore keys are in separate branches */ 1.1332 + return false; 1.1333 + } 1.1334 + } 1.1335 + } 1.1336 + /* return true because keys are in the same branch */ 1.1337 + return true; 1.1338 +} 1.1339 + 1.1340 +/* combine two keys into new key which covers both original keys */ 1.1341 +/* (result stored in first argument) */ 1.1342 +static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) { 1.1343 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 1.1344 + /* determine smallest depth */ 1.1345 + int depth1 = PGL_KEY_NODEDEPTH(dst); 1.1346 + int depth2 = PGL_KEY_NODEDEPTH(src); 1.1347 + int depth = (depth1 < depth2) ? depth1 : depth2; 1.1348 + /* check if keys are area keys (assuming that both keys have same type) */ 1.1349 + if (PGL_KEY_IS_AREAKEY(dst)) { 1.1350 + pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */ 1.1351 + int j = 0; /* bit offset for logarithmic object size bits */ 1.1352 + int k = 0; /* bit offset for latitude and longitude */ 1.1353 + /* fetch logarithmic object size information */ 1.1354 + int objsize1 = PGL_KEY_OBJSIZE(dst); 1.1355 + int objsize2 = PGL_KEY_OBJSIZE(src); 1.1356 + /* handle special cases for empty objects (universal and empty keys) */ 1.1357 + if ( 1.1358 + objsize1 > PGL_AREAKEY_MAXOBJSIZE || 1.1359 + objsize2 > PGL_AREAKEY_MAXOBJSIZE 1.1360 + ) { 1.1361 + if ( 1.1362 + objsize1 == PGL_KEY_OBJSIZE_EMPTY && 1.1363 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 1.1364 + ) pgl_key_set_empty(dst); 1.1365 + else pgl_key_set_universal(dst); 1.1366 + return; 1.1367 + } 1.1368 + /* iterate through key bits */ 1.1369 + for (i=0; i<depth; i++) { 1.1370 + /* every second bit is a bit describing the object size */ 1.1371 + if (i%2 == 0) { 1.1372 + /* increase bit counter for object size bits first */ 1.1373 + /* (handy when setting objsize variable) */ 1.1374 + j++; 1.1375 + /* check if object size bit is set in neither key */ 1.1376 + if (objsize1 >= j && objsize2 >= j) { 1.1377 + /* set objsize in destination buffer to indicate that size bit is 1.1378 + unset in destination buffer at the current bit position */ 1.1379 + dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j; 1.1380 + } 1.1381 + /* break if object size bit is set in one key only */ 1.1382 + else if (objsize1 >= j || objsize2 >= j) break; 1.1383 + } 1.1384 + /* all other bits describe latitude and longitude */ 1.1385 + else { 1.1386 + /* break if bit differs in both keys */ 1.1387 + if (PGL_KEY_LATLONBIT(dst, k)) { 1.1388 + if (!PGL_KEY_LATLONBIT(src, k)) break; 1.1389 + /* but set bit in destination buffer if bit is set in both keys */ 1.1390 + dstbuf[k/8] |= 0x80 >> (k%8); 1.1391 + } else if (PGL_KEY_LATLONBIT(src, k)) break; 1.1392 + /* increase bit counter for latitude/longitude bits */ 1.1393 + k++; 1.1394 + } 1.1395 + } 1.1396 + /* set common node depth and type bit (type bit = 1) */ 1.1397 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i; 1.1398 + /* copy contents of destination buffer to first key */ 1.1399 + memcpy(dst, dstbuf, sizeof(pgl_areakey)); 1.1400 + } 1.1401 + /* if not, keys are point keys */ 1.1402 + else { 1.1403 + pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */ 1.1404 + /* iterate through key bits */ 1.1405 + for (i=0; i<depth; i++) { 1.1406 + /* break if bit differs in both keys */ 1.1407 + if (PGL_KEY_LATLONBIT(dst, i)) { 1.1408 + if (!PGL_KEY_LATLONBIT(src, i)) break; 1.1409 + /* but set bit in destination buffer if bit is set in both keys */ 1.1410 + dstbuf[i/8] |= 0x80 >> (i%8); 1.1411 + } else if (PGL_KEY_LATLONBIT(src, i)) break; 1.1412 + } 1.1413 + /* set common node depth (type bit = 0) */ 1.1414 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i; 1.1415 + /* copy contents of destination buffer to first key */ 1.1416 + memcpy(dst, dstbuf, sizeof(pgl_pointkey)); 1.1417 + } 1.1418 +} 1.1419 + 1.1420 +/* determine center(!) boundaries and radius estimation of index key */ 1.1421 +static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) { 1.1422 + int i; 1.1423 + /* determine node depth */ 1.1424 + int depth = PGL_KEY_NODEDEPTH(key); 1.1425 + /* center point of possible result */ 1.1426 + double lat = 0; 1.1427 + double lon = 0; 1.1428 + /* maximum distance of real center point from key center */ 1.1429 + double dlat = 90; 1.1430 + double dlon = 180; 1.1431 + /* maximum radius of contained objects */ 1.1432 + double radius = 0; /* always return zero for point index keys */ 1.1433 + /* check if key is area key */ 1.1434 + if (PGL_KEY_IS_AREAKEY(key)) { 1.1435 + /* get logarithmic object size */ 1.1436 + int objsize = PGL_KEY_OBJSIZE(key); 1.1437 + /* handle special cases for empty objects (universal and empty keys) */ 1.1438 + if (objsize == PGL_KEY_OBJSIZE_EMPTY) { 1.1439 + pgl_box_set_empty(box); 1.1440 + return 0; 1.1441 + } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) { 1.1442 + box->lat_min = -90; 1.1443 + box->lat_max = 90; 1.1444 + box->lon_min = -180; 1.1445 + box->lon_max = 180; 1.1446 + return 0; /* any value >= 0 would do */ 1.1447 + } 1.1448 + /* calculate maximum possible radius of objects covered by the given key */ 1.1449 + if (objsize == 0) radius = INFINITY; 1.1450 + else { 1.1451 + radius = PGL_AREAKEY_REFOBJSIZE; 1.1452 + while (--objsize) radius /= M_SQRT2; 1.1453 + } 1.1454 + /* iterate over latitude and longitude bits in key */ 1.1455 + /* (every second bit is a latitude or longitude bit) */ 1.1456 + for (i=0; i<depth/2; i++) { 1.1457 + /* check if latitude bit */ 1.1458 + if (i%2 == 0) { 1.1459 + /* cut latitude dimension in half */ 1.1460 + dlat /= 2; 1.1461 + /* increase center latitude if bit is 1, otherwise decrease */ 1.1462 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 1.1463 + else lat -= dlat; 1.1464 + } 1.1465 + /* otherwise longitude bit */ 1.1466 + else { 1.1467 + /* cut longitude dimension in half */ 1.1468 + dlon /= 2; 1.1469 + /* increase center longitude if bit is 1, otherwise decrease */ 1.1470 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 1.1471 + else lon -= dlon; 1.1472 + } 1.1473 + } 1.1474 + } 1.1475 + /* if not, keys are point keys */ 1.1476 + else { 1.1477 + /* iterate over all bits in key */ 1.1478 + for (i=0; i<depth; i++) { 1.1479 + /* check if latitude bit */ 1.1480 + if (i%2 == 0) { 1.1481 + /* cut latitude dimension in half */ 1.1482 + dlat /= 2; 1.1483 + /* increase center latitude if bit is 1, otherwise decrease */ 1.1484 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 1.1485 + else lat -= dlat; 1.1486 + } 1.1487 + /* otherwise longitude bit */ 1.1488 + else { 1.1489 + /* cut longitude dimension in half */ 1.1490 + dlon /= 2; 1.1491 + /* increase center longitude if bit is 1, otherwise decrease */ 1.1492 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 1.1493 + else lon -= dlon; 1.1494 + } 1.1495 + } 1.1496 + } 1.1497 + /* calculate boundaries from center point and remaining dlat and dlon */ 1.1498 + /* (return values through pointer to box) */ 1.1499 + box->lat_min = lat - dlat; 1.1500 + box->lat_max = lat + dlat; 1.1501 + box->lon_min = lon - dlon; 1.1502 + box->lon_max = lon + dlon; 1.1503 + /* return radius (as a function return value) */ 1.1504 + return radius; 1.1505 +} 1.1506 + 1.1507 +/* estimator function for distance between point and index key */ 1.1508 +/* always returns a smaller value than actually correct or zero */ 1.1509 +static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) { 1.1510 + pgl_box box; /* center(!) bounding box of area index key */ 1.1511 + /* calculate center(!) bounding box and maximum radius of objects covered 1.1512 + by area index key (radius is zero for point index keys) */ 1.1513 + double distance = pgl_key_to_box(key, &box); 1.1514 + /* calculate estimated distance between bounding box of center point of 1.1515 + indexed object and point passed as second argument, then substract maximum 1.1516 + radius of objects covered by index key */ 1.1517 + distance = pgl_estimate_point_box_distance(point, &box) - distance; 1.1518 + /* truncate negative results to zero */ 1.1519 + if (distance <= 0) distance = 0; 1.1520 + /* return result */ 1.1521 + return distance; 1.1522 +} 1.1523 + 1.1524 + 1.1525 +/*---------------------------------* 1.1526 + * helper functions for text I/O * 1.1527 + *---------------------------------*/ 1.1528 + 1.1529 +#define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */ 1.1530 + 1.1531 +/* convert floating point number to string (round-trip safe) */ 1.1532 +static void pgl_print_float(char *buf, double flt) { 1.1533 + /* check if number is integral */ 1.1534 + if (trunc(flt) == flt) { 1.1535 + /* for integral floats use maximum precision */ 1.1536 + snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 1.1537 + } else { 1.1538 + /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */ 1.1539 + snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt); 1.1540 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt); 1.1541 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 1.1542 + } 1.1543 +} 1.1544 + 1.1545 +/* convert latitude floating point number (in degrees) to string */ 1.1546 +static void pgl_print_lat(char *buf, double lat) { 1.1547 + if (signbit(lat)) { 1.1548 + /* treat negative latitudes (including -0) as south */ 1.1549 + snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat); 1.1550 + } else { 1.1551 + /* treat positive latitudes (including +0) as north */ 1.1552 + snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat); 1.1553 + } 1.1554 +} 1.1555 + 1.1556 +/* convert longitude floating point number (in degrees) to string */ 1.1557 +static void pgl_print_lon(char *buf, double lon) { 1.1558 + if (signbit(lon)) { 1.1559 + /* treat negative longitudes (including -0) as west */ 1.1560 + snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon); 1.1561 + } else { 1.1562 + /* treat positive longitudes (including +0) as east */ 1.1563 + snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon); 1.1564 + } 1.1565 +} 1.1566 + 1.1567 +/* bit masks used as return value of pgl_scan() function */ 1.1568 +#define PGL_SCAN_NONE 0 /* no value has been parsed */ 1.1569 +#define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */ 1.1570 +#define PGL_SCAN_LON (1<<1) /* longitude has been parsed */ 1.1571 +#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */ 1.1572 + 1.1573 +/* parse a coordinate (can be latitude or longitude) */ 1.1574 +static int pgl_scan(char **str, double *lat, double *lon) { 1.1575 + double val; 1.1576 + int len; 1.1577 + if ( 1.1578 + sscanf(*str, " N %lf %n", &val, &len) || 1.1579 + sscanf(*str, " n %lf %n", &val, &len) 1.1580 + ) { 1.1581 + *str += len; *lat = val; return PGL_SCAN_LAT; 1.1582 + } 1.1583 + if ( 1.1584 + sscanf(*str, " S %lf %n", &val, &len) || 1.1585 + sscanf(*str, " s %lf %n", &val, &len) 1.1586 + ) { 1.1587 + *str += len; *lat = -val; return PGL_SCAN_LAT; 1.1588 + } 1.1589 + if ( 1.1590 + sscanf(*str, " E %lf %n", &val, &len) || 1.1591 + sscanf(*str, " e %lf %n", &val, &len) 1.1592 + ) { 1.1593 + *str += len; *lon = val; return PGL_SCAN_LON; 1.1594 + } 1.1595 + if ( 1.1596 + sscanf(*str, " W %lf %n", &val, &len) || 1.1597 + sscanf(*str, " w %lf %n", &val, &len) 1.1598 + ) { 1.1599 + *str += len; *lon = -val; return PGL_SCAN_LON; 1.1600 + } 1.1601 + return PGL_SCAN_NONE; 1.1602 +} 1.1603 + 1.1604 + 1.1605 +/*-----------------* 1.1606 + * SQL functions * 1.1607 + *-----------------*/ 1.1608 + 1.1609 +/* Note: These function names use "epoint", "ebox", etc. notation here instead 1.1610 + of "point", "box", etc. in order to distinguish them from any previously 1.1611 + defined functions. */ 1.1612 + 1.1613 +/* function needed for dummy types and/or not implemented features */ 1.1614 +PG_FUNCTION_INFO_V1(pgl_notimpl); 1.1615 +Datum pgl_notimpl(PG_FUNCTION_ARGS) { 1.1616 + ereport(ERROR, (errmsg("not implemented by pgLatLon"))); 1.1617 +} 1.1618 + 1.1619 +/* set point to latitude and longitude (including checks) */ 1.1620 +static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) { 1.1621 + /* reject infinite or NaN values */ 1.1622 + if (!isfinite(lat) || !isfinite(lon)) { 1.1623 + ereport(ERROR, ( 1.1624 + errcode(ERRCODE_DATA_EXCEPTION), 1.1625 + errmsg("epoint requires finite coordinates") 1.1626 + )); 1.1627 + } 1.1628 + /* check latitude bounds */ 1.1629 + if (lat < -90) { 1.1630 + ereport(WARNING, (errmsg("latitude exceeds south pole"))); 1.1631 + lat = -90; 1.1632 + } else if (lat > 90) { 1.1633 + ereport(WARNING, (errmsg("latitude exceeds north pole"))); 1.1634 + lat = 90; 1.1635 + } 1.1636 + /* check longitude bounds */ 1.1637 + if (lon < -180) { 1.1638 + ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized"))); 1.1639 + lon += 360 - trunc(lon / 360) * 360; 1.1640 + } else if (lon > 180) { 1.1641 + ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized"))); 1.1642 + lon -= 360 + trunc(lon / 360) * 360; 1.1643 + } 1.1644 + /* store rounded latitude/longitude values for round-trip safety */ 1.1645 + point->lat = pgl_round(lat); 1.1646 + point->lon = pgl_round(lon); 1.1647 +} 1.1648 + 1.1649 +/* create point ("epoint" in SQL) from latitude and longitude */ 1.1650 +PG_FUNCTION_INFO_V1(pgl_create_epoint); 1.1651 +Datum pgl_create_epoint(PG_FUNCTION_ARGS) { 1.1652 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 1.1653 + pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1)); 1.1654 + PG_RETURN_POINTER(point); 1.1655 +} 1.1656 + 1.1657 +/* parse point ("epoint" in SQL) */ 1.1658 +/* format: '[NS]<float> [EW]<float>' */ 1.1659 +PG_FUNCTION_INFO_V1(pgl_epoint_in); 1.1660 +Datum pgl_epoint_in(PG_FUNCTION_ARGS) { 1.1661 + char *str = PG_GETARG_CSTRING(0); /* input string */ 1.1662 + char *strptr = str; /* current position within string */ 1.1663 + int done = 0; /* bit mask storing if latitude or longitude was read */ 1.1664 + double lat, lon; /* parsed values as double precision floats */ 1.1665 + pgl_point *point; /* return value (to be palloc'ed) */ 1.1666 + /* parse two floats (each latitude or longitude) separated by white-space */ 1.1667 + done |= pgl_scan(&strptr, &lat, &lon); 1.1668 + if (strptr != str && isspace(strptr[-1])) { 1.1669 + done |= pgl_scan(&strptr, &lat, &lon); 1.1670 + } 1.1671 + /* require end of string, and latitude and longitude parsed successfully */ 1.1672 + if (strptr[0] || done != PGL_SCAN_LATLON) { 1.1673 + ereport(ERROR, ( 1.1674 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.1675 + errmsg("invalid input syntax for type epoint: \"%s\"", str) 1.1676 + )); 1.1677 + } 1.1678 + /* allocate memory for result */ 1.1679 + point = (pgl_point *)palloc(sizeof(pgl_point)); 1.1680 + /* set latitude and longitude (and perform checks) */ 1.1681 + pgl_epoint_set_latlon(point, lat, lon); 1.1682 + /* return result */ 1.1683 + PG_RETURN_POINTER(point); 1.1684 +} 1.1685 + 1.1686 +/* set sample count for numerical integration (including checks) */ 1.1687 +static void pgl_epoint_set_sample_count(pgl_point_sc *search, int32 samples) { 1.1688 + /* require minimum of 6 samples */ 1.1689 + if (samples < 6) { 1.1690 + ereport(ERROR, ( 1.1691 + errcode(ERRCODE_DATA_EXCEPTION), 1.1692 + errmsg("too few sample points for numerical integration (minimum 6)") 1.1693 + )); 1.1694 + } 1.1695 + /* limit sample count to avoid integer overflows on memory allocation */ 1.1696 + if (samples > PGL_CLUSTER_MAXPOINTS) { 1.1697 + ereport(ERROR, ( 1.1698 + errcode(ERRCODE_DATA_EXCEPTION), 1.1699 + errmsg( 1.1700 + "too many sample points for numerical integration (maximum %i)", 1.1701 + PGL_CLUSTER_MAXPOINTS 1.1702 + ) 1.1703 + )); 1.1704 + } 1.1705 + search->samples = samples; 1.1706 +} 1.1707 + 1.1708 +/* create point with sample count for fair distance calculation 1.1709 + ("epoint_with_sample_count" in SQL) from epoint and integer */ 1.1710 +PG_FUNCTION_INFO_V1(pgl_create_epoint_with_sample_count); 1.1711 +Datum pgl_create_epoint_with_sample_count(PG_FUNCTION_ARGS) { 1.1712 + pgl_point_sc *search = (pgl_point_sc *)palloc(sizeof(pgl_point_sc)); 1.1713 + search->point = *(pgl_point *)PG_GETARG_POINTER(0); 1.1714 + pgl_epoint_set_sample_count(search, PG_GETARG_INT32(1)); 1.1715 + PG_RETURN_POINTER(search); 1.1716 +} 1.1717 + 1.1718 +/* parse point with sample count ("epoint_with_sample_count" in SQL) */ 1.1719 +/* format: '[NS]<float> [EW]<float> <integer>' */ 1.1720 +PG_FUNCTION_INFO_V1(pgl_epoint_with_sample_count_in); 1.1721 +Datum pgl_epoint_with_sample_count_in(PG_FUNCTION_ARGS) { 1.1722 + char *str = PG_GETARG_CSTRING(0); /* input string */ 1.1723 + char *strptr = str; /* current position within string */ 1.1724 + double lat, lon; /* parsed values for latitude and longitude */ 1.1725 + int samples; /* parsed value for sample count */ 1.1726 + int valid = 0; /* number of valid chars */ 1.1727 + int done = 0; /* stores if latitude and/or longitude was read */ 1.1728 + pgl_point_sc *search; /* return value (to be palloc'ed) */ 1.1729 + /* demand three blocks separated by whitespace */ 1.1730 + sscanf(strptr, " %*s %*s %*s %n", &valid); 1.1731 + /* if three blocks separated by whitespace exist, parse those blocks */ 1.1732 + if (strptr[valid] == 0) { 1.1733 + /* parse latitude and longitude */ 1.1734 + done |= pgl_scan(&strptr, &lat, &lon); 1.1735 + done |= pgl_scan(&strptr, &lat, &lon); 1.1736 + /* parse sample count (while incr. strptr by number of bytes parsed) */ 1.1737 + valid = 0; 1.1738 + if (sscanf(strptr, " %d %n", &samples, &valid) == 1) strptr += valid; 1.1739 + } 1.1740 + /* require end of string and both latitude and longitude being parsed */ 1.1741 + if (strptr[0] || done != PGL_SCAN_LATLON) { 1.1742 + ereport(ERROR, ( 1.1743 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.1744 + errmsg("invalid input syntax for type ecircle: \"%s\"", str) 1.1745 + )); 1.1746 + } 1.1747 + /* allocate memory for result */ 1.1748 + search = (pgl_point_sc *)palloc(sizeof(pgl_point_sc)); 1.1749 + /* set latitude, longitude, and sample count (while performing checks) */ 1.1750 + pgl_epoint_set_latlon(&search->point, lat, lon); 1.1751 + pgl_epoint_set_sample_count(search, samples); 1.1752 + /* return result */ 1.1753 + PG_RETURN_POINTER(search); 1.1754 +} 1.1755 + 1.1756 +/* create box ("ebox" in SQL) that is empty */ 1.1757 +PG_FUNCTION_INFO_V1(pgl_create_empty_ebox); 1.1758 +Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) { 1.1759 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 1.1760 + pgl_box_set_empty(box); 1.1761 + PG_RETURN_POINTER(box); 1.1762 +} 1.1763 + 1.1764 +/* set box to given boundaries (including checks) */ 1.1765 +static void pgl_ebox_set_boundaries( 1.1766 + pgl_box *box, 1.1767 + double lat_min, double lat_max, double lon_min, double lon_max 1.1768 +) { 1.1769 + /* if minimum latitude is greater than maximum latitude, return empty box */ 1.1770 + if (lat_min > lat_max) { 1.1771 + pgl_box_set_empty(box); 1.1772 + return; 1.1773 + } 1.1774 + /* otherwise reject infinite or NaN values */ 1.1775 + if ( 1.1776 + !isfinite(lat_min) || !isfinite(lat_max) || 1.1777 + !isfinite(lon_min) || !isfinite(lon_max) 1.1778 + ) { 1.1779 + ereport(ERROR, ( 1.1780 + errcode(ERRCODE_DATA_EXCEPTION), 1.1781 + errmsg("ebox requires finite coordinates") 1.1782 + )); 1.1783 + } 1.1784 + /* check latitude bounds */ 1.1785 + if (lat_max < -90) { 1.1786 + ereport(WARNING, (errmsg("northern latitude exceeds south pole"))); 1.1787 + lat_max = -90; 1.1788 + } else if (lat_max > 90) { 1.1789 + ereport(WARNING, (errmsg("northern latitude exceeds north pole"))); 1.1790 + lat_max = 90; 1.1791 + } 1.1792 + if (lat_min < -90) { 1.1793 + ereport(WARNING, (errmsg("southern latitude exceeds south pole"))); 1.1794 + lat_min = -90; 1.1795 + } else if (lat_min > 90) { 1.1796 + ereport(WARNING, (errmsg("southern latitude exceeds north pole"))); 1.1797 + lat_min = 90; 1.1798 + } 1.1799 + /* check if all longitudes are included */ 1.1800 + if (lon_max - lon_min >= 360) { 1.1801 + if (lon_max - lon_min > 360) ereport(WARNING, ( 1.1802 + errmsg("longitude coverage greater than 360 degrees") 1.1803 + )); 1.1804 + lon_min = -180; 1.1805 + lon_max = 180; 1.1806 + } else { 1.1807 + /* normalize longitude bounds */ 1.1808 + if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360; 1.1809 + else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360; 1.1810 + if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360; 1.1811 + else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360; 1.1812 + } 1.1813 + /* store rounded latitude/longitude values for round-trip safety */ 1.1814 + box->lat_min = pgl_round(lat_min); 1.1815 + box->lat_max = pgl_round(lat_max); 1.1816 + box->lon_min = pgl_round(lon_min); 1.1817 + box->lon_max = pgl_round(lon_max); 1.1818 + /* ensure that rounding does not change orientation */ 1.1819 + if (lon_min > lon_max && box->lon_min == box->lon_max) { 1.1820 + box->lon_min = -180; 1.1821 + box->lon_max = 180; 1.1822 + } 1.1823 +} 1.1824 + 1.1825 +/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */ 1.1826 +PG_FUNCTION_INFO_V1(pgl_create_ebox); 1.1827 +Datum pgl_create_ebox(PG_FUNCTION_ARGS) { 1.1828 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 1.1829 + pgl_ebox_set_boundaries( 1.1830 + box, 1.1831 + PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), 1.1832 + PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3) 1.1833 + ); 1.1834 + PG_RETURN_POINTER(box); 1.1835 +} 1.1836 + 1.1837 +/* create box ("ebox" in SQL) from two points ("epoint"s) */ 1.1838 +/* (can not be used to cover a longitude range of more than 120 degrees) */ 1.1839 +PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints); 1.1840 +Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) { 1.1841 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 1.1842 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 1.1843 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 1.1844 + double lat_min, lat_max, lon_min, lon_max; 1.1845 + double dlon; /* longitude range (delta longitude) */ 1.1846 + /* order latitude and longitude boundaries */ 1.1847 + if (point2->lat < point1->lat) { 1.1848 + lat_min = point2->lat; 1.1849 + lat_max = point1->lat; 1.1850 + } else { 1.1851 + lat_min = point1->lat; 1.1852 + lat_max = point2->lat; 1.1853 + } 1.1854 + if (point2->lon < point1->lon) { 1.1855 + lon_min = point2->lon; 1.1856 + lon_max = point1->lon; 1.1857 + } else { 1.1858 + lon_min = point1->lon; 1.1859 + lon_max = point2->lon; 1.1860 + } 1.1861 + /* calculate longitude range (round to avoid floating point errors) */ 1.1862 + dlon = pgl_round(lon_max - lon_min); 1.1863 + /* determine east-west direction */ 1.1864 + if (dlon >= 240) { 1.1865 + /* assume that 180th meridian is crossed and swap min/max longitude */ 1.1866 + double swap = lon_min; lon_min = lon_max; lon_max = swap; 1.1867 + } else if (dlon > 120) { 1.1868 + /* unclear orientation since delta longitude > 120 */ 1.1869 + ereport(ERROR, ( 1.1870 + errcode(ERRCODE_DATA_EXCEPTION), 1.1871 + errmsg("can not determine east/west orientation for ebox") 1.1872 + )); 1.1873 + } 1.1874 + /* use boundaries to setup box (and perform checks) */ 1.1875 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 1.1876 + /* return result */ 1.1877 + PG_RETURN_POINTER(box); 1.1878 +} 1.1879 + 1.1880 +/* parse box ("ebox" in SQL) */ 1.1881 +/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>' 1.1882 + or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */ 1.1883 +PG_FUNCTION_INFO_V1(pgl_ebox_in); 1.1884 +Datum pgl_ebox_in(PG_FUNCTION_ARGS) { 1.1885 + char *str = PG_GETARG_CSTRING(0); /* input string */ 1.1886 + char *str_lower; /* lower case version of input string */ 1.1887 + char *strptr; /* current position within string */ 1.1888 + int valid; /* number of valid chars */ 1.1889 + int done; /* specifies if latitude or longitude was read */ 1.1890 + double val; /* temporary variable */ 1.1891 + int lat_count = 0; /* count of latitude values parsed */ 1.1892 + int lon_count = 0; /* count of longitufde values parsed */ 1.1893 + double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */ 1.1894 + pgl_box *box; /* return value (to be palloc'ed) */ 1.1895 + /* lowercase input */ 1.1896 + str_lower = psprintf("%s", str); 1.1897 + for (strptr=str_lower; *strptr; strptr++) { 1.1898 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 1.1899 + } 1.1900 + /* reset reading position to start of (lowercase) string */ 1.1901 + strptr = str_lower; 1.1902 + /* check if empty box */ 1.1903 + valid = 0; 1.1904 + sscanf(strptr, " empty %n", &valid); 1.1905 + if (valid && strptr[valid] == 0) { 1.1906 + /* allocate and return empty box */ 1.1907 + box = (pgl_box *)palloc(sizeof(pgl_box)); 1.1908 + pgl_box_set_empty(box); 1.1909 + PG_RETURN_POINTER(box); 1.1910 + } 1.1911 + /* demand four blocks separated by whitespace */ 1.1912 + valid = 0; 1.1913 + sscanf(strptr, " %*s %*s %*s %*s %n", &valid); 1.1914 + /* if four blocks separated by whitespace exist, parse those blocks */ 1.1915 + if (strptr[valid] == 0) while (strptr[0]) { 1.1916 + /* parse either latitude or longitude (whichever found in input string) */ 1.1917 + done = pgl_scan(&strptr, &val, &val); 1.1918 + /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */ 1.1919 + if (done == PGL_SCAN_LAT) { 1.1920 + if (!lat_count) lat_min = val; else lat_max = val; 1.1921 + lat_count++; 1.1922 + } else if (done == PGL_SCAN_LON) { 1.1923 + if (!lon_count) lon_min = val; else lon_max = val; 1.1924 + lon_count++; 1.1925 + } else { 1.1926 + break; 1.1927 + } 1.1928 + } 1.1929 + /* require end of string, and two latitude and two longitude values */ 1.1930 + if (strptr[0] || lat_count != 2 || lon_count != 2) { 1.1931 + ereport(ERROR, ( 1.1932 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.1933 + errmsg("invalid input syntax for type ebox: \"%s\"", str) 1.1934 + )); 1.1935 + } 1.1936 + /* free lower case string */ 1.1937 + pfree(str_lower); 1.1938 + /* order boundaries (maximum greater than minimum) */ 1.1939 + if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; } 1.1940 + if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; } 1.1941 + /* allocate memory for result */ 1.1942 + box = (pgl_box *)palloc(sizeof(pgl_box)); 1.1943 + /* set boundaries (and perform checks) */ 1.1944 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 1.1945 + /* return result */ 1.1946 + PG_RETURN_POINTER(box); 1.1947 +} 1.1948 + 1.1949 +/* set circle to given latitude, longitude, and radius (including checks) */ 1.1950 +static void pgl_ecircle_set_latlon_radius( 1.1951 + pgl_circle *circle, double lat, double lon, double radius 1.1952 +) { 1.1953 + /* set center point (including checks) */ 1.1954 + pgl_epoint_set_latlon(&(circle->center), lat, lon); 1.1955 + /* handle non-positive radius */ 1.1956 + if (isnan(radius)) { 1.1957 + ereport(ERROR, ( 1.1958 + errcode(ERRCODE_DATA_EXCEPTION), 1.1959 + errmsg("invalid radius for ecircle") 1.1960 + )); 1.1961 + } 1.1962 + if (radius == 0) radius = 0; /* avoids -0 */ 1.1963 + else if (radius < 0) { 1.1964 + if (isfinite(radius)) { 1.1965 + ereport(NOTICE, (errmsg("negative radius converted to minus infinity"))); 1.1966 + } 1.1967 + radius = -INFINITY; 1.1968 + } 1.1969 + /* store radius (round-trip safety is ensured by pgl_print_float) */ 1.1970 + circle->radius = radius; 1.1971 +} 1.1972 + 1.1973 +/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */ 1.1974 +PG_FUNCTION_INFO_V1(pgl_create_ecircle); 1.1975 +Datum pgl_create_ecircle(PG_FUNCTION_ARGS) { 1.1976 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 1.1977 + pgl_ecircle_set_latlon_radius( 1.1978 + circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2) 1.1979 + ); 1.1980 + PG_RETURN_POINTER(circle); 1.1981 +} 1.1982 + 1.1983 +/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */ 1.1984 +PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint); 1.1985 +Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) { 1.1986 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.1987 + double radius = PG_GETARG_FLOAT8(1); 1.1988 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 1.1989 + /* set latitude, longitude, radius (and perform checks) */ 1.1990 + pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius); 1.1991 + /* return result */ 1.1992 + PG_RETURN_POINTER(circle); 1.1993 +} 1.1994 + 1.1995 +/* parse circle ("ecircle" in SQL) */ 1.1996 +/* format: '[NS]<float> [EW]<float> <float>' */ 1.1997 +PG_FUNCTION_INFO_V1(pgl_ecircle_in); 1.1998 +Datum pgl_ecircle_in(PG_FUNCTION_ARGS) { 1.1999 + char *str = PG_GETARG_CSTRING(0); /* input string */ 1.2000 + char *strptr = str; /* current position within string */ 1.2001 + double lat, lon, radius; /* parsed values as double precision floats */ 1.2002 + int valid = 0; /* number of valid chars */ 1.2003 + int done = 0; /* stores if latitude and/or longitude was read */ 1.2004 + pgl_circle *circle; /* return value (to be palloc'ed) */ 1.2005 + /* demand three blocks separated by whitespace */ 1.2006 + sscanf(strptr, " %*s %*s %*s %n", &valid); 1.2007 + /* if three blocks separated by whitespace exist, parse those blocks */ 1.2008 + if (strptr[valid] == 0) { 1.2009 + /* parse latitude and longitude */ 1.2010 + done |= pgl_scan(&strptr, &lat, &lon); 1.2011 + done |= pgl_scan(&strptr, &lat, &lon); 1.2012 + /* parse radius (while incrementing strptr by number of bytes parsed) */ 1.2013 + valid = 0; 1.2014 + if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid; 1.2015 + } 1.2016 + /* require end of string and both latitude and longitude being parsed */ 1.2017 + if (strptr[0] || done != PGL_SCAN_LATLON) { 1.2018 + ereport(ERROR, ( 1.2019 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.2020 + errmsg("invalid input syntax for type ecircle: \"%s\"", str) 1.2021 + )); 1.2022 + } 1.2023 + /* allocate memory for result */ 1.2024 + circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 1.2025 + /* set latitude, longitude, radius (and perform checks) */ 1.2026 + pgl_ecircle_set_latlon_radius(circle, lat, lon, radius); 1.2027 + /* return result */ 1.2028 + PG_RETURN_POINTER(circle); 1.2029 +} 1.2030 + 1.2031 +/* parse cluster ("ecluster" in SQL) */ 1.2032 +PG_FUNCTION_INFO_V1(pgl_ecluster_in); 1.2033 +Datum pgl_ecluster_in(PG_FUNCTION_ARGS) { 1.2034 + int i; 1.2035 + char *str = PG_GETARG_CSTRING(0); /* input string */ 1.2036 + char *str_lower; /* lower case version of input string */ 1.2037 + char *strptr; /* pointer to current reading position of input */ 1.2038 + int npoints_total = 0; /* total number of points in cluster */ 1.2039 + int nentries = 0; /* total number of entries */ 1.2040 + pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */ 1.2041 + int entries_buflen = 4; /* maximum number of elements in entries array */ 1.2042 + int valid; /* number of valid chars processed */ 1.2043 + double lat, lon; /* latitude and longitude of parsed point */ 1.2044 + int entrytype; /* current entry type */ 1.2045 + int npoints; /* number of points in current entry */ 1.2046 + pgl_point *points; /* array of pgl_point for pgl_newentry */ 1.2047 + int points_buflen; /* maximum number of elements in points array */ 1.2048 + int done; /* return value of pgl_scan function */ 1.2049 + pgl_cluster *cluster; /* created cluster */ 1.2050 + /* lowercase input */ 1.2051 + str_lower = psprintf("%s", str); 1.2052 + for (strptr=str_lower; *strptr; strptr++) { 1.2053 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 1.2054 + } 1.2055 + /* reset reading position to start of (lowercase) string */ 1.2056 + strptr = str_lower; 1.2057 + /* allocate initial buffer for entries */ 1.2058 + entries = palloc(entries_buflen * sizeof(pgl_newentry)); 1.2059 + /* parse until end of string */ 1.2060 + while (strptr[0]) { 1.2061 + /* require previous white-space or closing parenthesis before next token */ 1.2062 + if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') { 1.2063 + goto pgl_ecluster_in_error; 1.2064 + } 1.2065 + /* ignore token "empty" */ 1.2066 + valid = 0; sscanf(strptr, " empty %n", &valid); 1.2067 + if (valid) { strptr += valid; continue; } 1.2068 + /* test for "point" token */ 1.2069 + valid = 0; sscanf(strptr, " point ( %n", &valid); 1.2070 + if (valid) { 1.2071 + strptr += valid; 1.2072 + entrytype = PGL_ENTRY_POINT; 1.2073 + goto pgl_ecluster_in_type_ok; 1.2074 + } 1.2075 + /* test for "path" token */ 1.2076 + valid = 0; sscanf(strptr, " path ( %n", &valid); 1.2077 + if (valid) { 1.2078 + strptr += valid; 1.2079 + entrytype = PGL_ENTRY_PATH; 1.2080 + goto pgl_ecluster_in_type_ok; 1.2081 + } 1.2082 + /* test for "outline" token */ 1.2083 + valid = 0; sscanf(strptr, " outline ( %n", &valid); 1.2084 + if (valid) { 1.2085 + strptr += valid; 1.2086 + entrytype = PGL_ENTRY_OUTLINE; 1.2087 + goto pgl_ecluster_in_type_ok; 1.2088 + } 1.2089 + /* test for "polygon" token */ 1.2090 + valid = 0; sscanf(strptr, " polygon ( %n", &valid); 1.2091 + if (valid) { 1.2092 + strptr += valid; 1.2093 + entrytype = PGL_ENTRY_POLYGON; 1.2094 + goto pgl_ecluster_in_type_ok; 1.2095 + } 1.2096 + /* error if no valid token found */ 1.2097 + goto pgl_ecluster_in_error; 1.2098 + pgl_ecluster_in_type_ok: 1.2099 + /* check if pgl_newentry array needs to grow */ 1.2100 + if (nentries == entries_buflen) { 1.2101 + pgl_newentry *newbuf; 1.2102 + entries_buflen *= 2; 1.2103 + newbuf = palloc(entries_buflen * sizeof(pgl_newentry)); 1.2104 + memcpy(newbuf, entries, nentries * sizeof(pgl_newentry)); 1.2105 + pfree(entries); 1.2106 + entries = newbuf; 1.2107 + } 1.2108 + /* reset number of points for current entry */ 1.2109 + npoints = 0; 1.2110 + /* allocate array for points */ 1.2111 + points_buflen = 4; 1.2112 + points = palloc(points_buflen * sizeof(pgl_point)); 1.2113 + /* parse until closing parenthesis */ 1.2114 + while (strptr[0] != ')') { 1.2115 + /* error on unexpected end of string */ 1.2116 + if (strptr[0] == 0) goto pgl_ecluster_in_error; 1.2117 + /* mark neither latitude nor longitude as read */ 1.2118 + done = PGL_SCAN_NONE; 1.2119 + /* require white-space before second, third, etc. point */ 1.2120 + if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 1.2121 + /* scan latitude (or longitude) */ 1.2122 + done |= pgl_scan(&strptr, &lat, &lon); 1.2123 + /* require white-space before second coordinate */ 1.2124 + if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 1.2125 + /* scan longitude (or latitude) */ 1.2126 + done |= pgl_scan(&strptr, &lat, &lon); 1.2127 + /* error unless both latitude and longitude were parsed */ 1.2128 + if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error; 1.2129 + /* throw error if number of points is too high */ 1.2130 + if (npoints_total == PGL_CLUSTER_MAXPOINTS) { 1.2131 + ereport(ERROR, ( 1.2132 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.2133 + errmsg( 1.2134 + "too many points for ecluster entry (maximum %i)", 1.2135 + PGL_CLUSTER_MAXPOINTS 1.2136 + ) 1.2137 + )); 1.2138 + } 1.2139 + /* check if pgl_point array needs to grow */ 1.2140 + if (npoints == points_buflen) { 1.2141 + pgl_point *newbuf; 1.2142 + points_buflen *= 2; 1.2143 + newbuf = palloc(points_buflen * sizeof(pgl_point)); 1.2144 + memcpy(newbuf, points, npoints * sizeof(pgl_point)); 1.2145 + pfree(points); 1.2146 + points = newbuf; 1.2147 + } 1.2148 + /* append point to pgl_point array (includes checks) */ 1.2149 + pgl_epoint_set_latlon(&(points[npoints++]), lat, lon); 1.2150 + /* increase total number of points */ 1.2151 + npoints_total++; 1.2152 + } 1.2153 + /* error if entry has no points */ 1.2154 + if (!npoints) goto pgl_ecluster_in_error; 1.2155 + /* entries with one point are automatically of type "point" */ 1.2156 + if (npoints == 1) entrytype = PGL_ENTRY_POINT; 1.2157 + /* if entries have more than one point */ 1.2158 + else { 1.2159 + /* throw error if entry type is "point" */ 1.2160 + if (entrytype == PGL_ENTRY_POINT) { 1.2161 + ereport(ERROR, ( 1.2162 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.2163 + errmsg("invalid input syntax for type ecluster (point entry with more than one point)") 1.2164 + )); 1.2165 + } 1.2166 + /* coerce outlines and polygons with more than 2 points to be a path */ 1.2167 + if (npoints == 2) entrytype = PGL_ENTRY_PATH; 1.2168 + } 1.2169 + /* append entry to pgl_newentry array */ 1.2170 + entries[nentries].entrytype = entrytype; 1.2171 + entries[nentries].npoints = npoints; 1.2172 + entries[nentries].points = points; 1.2173 + nentries++; 1.2174 + /* consume closing parenthesis */ 1.2175 + strptr++; 1.2176 + /* consume white-space */ 1.2177 + while (isspace(strptr[0])) strptr++; 1.2178 + } 1.2179 + /* free lower case string */ 1.2180 + pfree(str_lower); 1.2181 + /* create cluster from pgl_newentry array */ 1.2182 + cluster = pgl_new_cluster(nentries, entries); 1.2183 + /* free pgl_newentry array */ 1.2184 + for (i=0; i<nentries; i++) pfree(entries[i].points); 1.2185 + pfree(entries); 1.2186 + /* set bounding circle of cluster and check east/west orientation */ 1.2187 + if (!pgl_finalize_cluster(cluster)) { 1.2188 + ereport(ERROR, ( 1.2189 + errcode(ERRCODE_DATA_EXCEPTION), 1.2190 + errmsg("can not determine east/west orientation for ecluster"), 1.2191 + errhint("Ensure that each entry has a longitude span of less than 180 degrees.") 1.2192 + )); 1.2193 + } 1.2194 + /* return cluster */ 1.2195 + PG_RETURN_POINTER(cluster); 1.2196 + /* code to throw error */ 1.2197 + pgl_ecluster_in_error: 1.2198 + ereport(ERROR, ( 1.2199 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 1.2200 + errmsg("invalid input syntax for type ecluster: \"%s\"", str) 1.2201 + )); 1.2202 +} 1.2203 + 1.2204 +/* convert point ("epoint") to string representation */ 1.2205 +PG_FUNCTION_INFO_V1(pgl_epoint_out); 1.2206 +Datum pgl_epoint_out(PG_FUNCTION_ARGS) { 1.2207 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2208 + char latstr[PGL_NUMBUFLEN]; 1.2209 + char lonstr[PGL_NUMBUFLEN]; 1.2210 + pgl_print_lat(latstr, point->lat); 1.2211 + pgl_print_lon(lonstr, point->lon); 1.2212 + PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr)); 1.2213 +} 1.2214 + 1.2215 +/* convert point with sample count ("epoint_with_sample_count") to str. rep. */ 1.2216 +PG_FUNCTION_INFO_V1(pgl_epoint_with_sample_count_out); 1.2217 +Datum pgl_epoint_with_sample_count_out(PG_FUNCTION_ARGS) { 1.2218 + pgl_point_sc *search = (pgl_point_sc *)PG_GETARG_POINTER(0); 1.2219 + char latstr[PGL_NUMBUFLEN]; 1.2220 + char lonstr[PGL_NUMBUFLEN]; 1.2221 + pgl_print_lat(latstr, search->point.lat); 1.2222 + pgl_print_lon(lonstr, search->point.lon); 1.2223 + PG_RETURN_CSTRING( 1.2224 + psprintf("%s %s %i", latstr, lonstr, (int)search->samples) 1.2225 + ); 1.2226 +} 1.2227 + 1.2228 +/* convert box ("ebox") to string representation */ 1.2229 +PG_FUNCTION_INFO_V1(pgl_ebox_out); 1.2230 +Datum pgl_ebox_out(PG_FUNCTION_ARGS) { 1.2231 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 1.2232 + double lon_min = box->lon_min; 1.2233 + double lon_max = box->lon_max; 1.2234 + char lat_min_str[PGL_NUMBUFLEN]; 1.2235 + char lat_max_str[PGL_NUMBUFLEN]; 1.2236 + char lon_min_str[PGL_NUMBUFLEN]; 1.2237 + char lon_max_str[PGL_NUMBUFLEN]; 1.2238 + /* return string "empty" if box is set to be empty */ 1.2239 + if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty"); 1.2240 + /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */ 1.2241 + /* (required since pgl_box_in orders the longitude boundaries) */ 1.2242 + if (lon_min > lon_max) { 1.2243 + if (lon_min + lon_max >= 0) lon_min -= 360; 1.2244 + else lon_max += 360; 1.2245 + } 1.2246 + /* format and return result */ 1.2247 + pgl_print_lat(lat_min_str, box->lat_min); 1.2248 + pgl_print_lat(lat_max_str, box->lat_max); 1.2249 + pgl_print_lon(lon_min_str, lon_min); 1.2250 + pgl_print_lon(lon_max_str, lon_max); 1.2251 + PG_RETURN_CSTRING(psprintf( 1.2252 + "%s %s %s %s", 1.2253 + lat_min_str, lon_min_str, lat_max_str, lon_max_str 1.2254 + )); 1.2255 +} 1.2256 + 1.2257 +/* convert circle ("ecircle") to string representation */ 1.2258 +PG_FUNCTION_INFO_V1(pgl_ecircle_out); 1.2259 +Datum pgl_ecircle_out(PG_FUNCTION_ARGS) { 1.2260 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 1.2261 + char latstr[PGL_NUMBUFLEN]; 1.2262 + char lonstr[PGL_NUMBUFLEN]; 1.2263 + char radstr[PGL_NUMBUFLEN]; 1.2264 + pgl_print_lat(latstr, circle->center.lat); 1.2265 + pgl_print_lon(lonstr, circle->center.lon); 1.2266 + pgl_print_float(radstr, circle->radius); 1.2267 + PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr)); 1.2268 +} 1.2269 + 1.2270 +/* convert cluster ("ecluster") to string representation */ 1.2271 +PG_FUNCTION_INFO_V1(pgl_ecluster_out); 1.2272 +Datum pgl_ecluster_out(PG_FUNCTION_ARGS) { 1.2273 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2274 + char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */ 1.2275 + char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */ 1.2276 + char ***strings; /* array of array of strings */ 1.2277 + char *string; /* string of current token */ 1.2278 + char *res, *resptr; /* result and pointer to current write position */ 1.2279 + size_t reslen = 1; /* length of result (init with 1 for terminator) */ 1.2280 + int npoints; /* number of points of current entry */ 1.2281 + int i, j; /* i: entry, j: point in entry */ 1.2282 + /* handle empty clusters */ 1.2283 + if (cluster->nentries == 0) { 1.2284 + /* free detoasted cluster (if copy) */ 1.2285 + PG_FREE_IF_COPY(cluster, 0); 1.2286 + /* return static result */ 1.2287 + PG_RETURN_CSTRING("empty"); 1.2288 + } 1.2289 + /* allocate array of array of strings */ 1.2290 + strings = palloc(cluster->nentries * sizeof(char **)); 1.2291 + /* iterate over all entries in cluster */ 1.2292 + for (i=0; i<cluster->nentries; i++) { 1.2293 + /* get number of points in entry */ 1.2294 + npoints = cluster->entries[i].npoints; 1.2295 + /* allocate array of strings (one string for each point plus two extra) */ 1.2296 + strings[i] = palloc((2 + npoints) * sizeof(char *)); 1.2297 + /* determine opening string */ 1.2298 + switch (cluster->entries[i].entrytype) { 1.2299 + case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break; 1.2300 + case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break; 1.2301 + case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break; 1.2302 + case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break; 1.2303 + default: string = (i==0)?"unknown" :" unknown"; 1.2304 + } 1.2305 + /* use opening string as first string in array */ 1.2306 + strings[i][0] = string; 1.2307 + /* update result length (for allocating result string later) */ 1.2308 + reslen += strlen(string); 1.2309 + /* iterate over all points */ 1.2310 + for (j=0; j<npoints; j++) { 1.2311 + /* create string representation of point */ 1.2312 + pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat); 1.2313 + pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon); 1.2314 + string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr); 1.2315 + /* copy string pointer to string array */ 1.2316 + strings[i][j+1] = string; 1.2317 + /* update result length (for allocating result string later) */ 1.2318 + reslen += strlen(string); 1.2319 + } 1.2320 + /* use closing parenthesis as last string in array */ 1.2321 + strings[i][npoints+1] = ")"; 1.2322 + /* update result length (for allocating result string later) */ 1.2323 + reslen++; 1.2324 + } 1.2325 + /* allocate result string */ 1.2326 + res = palloc(reslen); 1.2327 + /* set write pointer to begin of result string */ 1.2328 + resptr = res; 1.2329 + /* copy strings into result string */ 1.2330 + for (i=0; i<cluster->nentries; i++) { 1.2331 + npoints = cluster->entries[i].npoints; 1.2332 + for (j=0; j<npoints+2; j++) { 1.2333 + string = strings[i][j]; 1.2334 + strcpy(resptr, string); 1.2335 + resptr += strlen(string); 1.2336 + /* free strings allocated by psprintf */ 1.2337 + if (j != 0 && j != npoints+1) pfree(string); 1.2338 + } 1.2339 + /* free array of strings */ 1.2340 + pfree(strings[i]); 1.2341 + } 1.2342 + /* free array of array of strings */ 1.2343 + pfree(strings); 1.2344 + /* free detoasted cluster (if copy) */ 1.2345 + PG_FREE_IF_COPY(cluster, 0); 1.2346 + /* return result */ 1.2347 + PG_RETURN_CSTRING(res); 1.2348 +} 1.2349 + 1.2350 +/* binary input function for point ("epoint") */ 1.2351 +PG_FUNCTION_INFO_V1(pgl_epoint_recv); 1.2352 +Datum pgl_epoint_recv(PG_FUNCTION_ARGS) { 1.2353 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 1.2354 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 1.2355 + point->lat = pq_getmsgfloat8(buf); 1.2356 + point->lon = pq_getmsgfloat8(buf); 1.2357 + PG_RETURN_POINTER(point); 1.2358 +} 1.2359 + 1.2360 +/* binary input function for box ("ebox") */ 1.2361 +PG_FUNCTION_INFO_V1(pgl_ebox_recv); 1.2362 +Datum pgl_ebox_recv(PG_FUNCTION_ARGS) { 1.2363 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 1.2364 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 1.2365 + box->lat_min = pq_getmsgfloat8(buf); 1.2366 + box->lat_max = pq_getmsgfloat8(buf); 1.2367 + box->lon_min = pq_getmsgfloat8(buf); 1.2368 + box->lon_max = pq_getmsgfloat8(buf); 1.2369 + PG_RETURN_POINTER(box); 1.2370 +} 1.2371 + 1.2372 +/* binary input function for circle ("ecircle") */ 1.2373 +PG_FUNCTION_INFO_V1(pgl_ecircle_recv); 1.2374 +Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) { 1.2375 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 1.2376 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 1.2377 + circle->center.lat = pq_getmsgfloat8(buf); 1.2378 + circle->center.lon = pq_getmsgfloat8(buf); 1.2379 + circle->radius = pq_getmsgfloat8(buf); 1.2380 + PG_RETURN_POINTER(circle); 1.2381 +} 1.2382 + 1.2383 +/* TODO: binary receive function for cluster */ 1.2384 + 1.2385 +/* binary output function for point ("epoint") */ 1.2386 +PG_FUNCTION_INFO_V1(pgl_epoint_send); 1.2387 +Datum pgl_epoint_send(PG_FUNCTION_ARGS) { 1.2388 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2389 + StringInfoData buf; 1.2390 + pq_begintypsend(&buf); 1.2391 + pq_sendfloat8(&buf, point->lat); 1.2392 + pq_sendfloat8(&buf, point->lon); 1.2393 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 1.2394 +} 1.2395 + 1.2396 +/* binary output function for box ("ebox") */ 1.2397 +PG_FUNCTION_INFO_V1(pgl_ebox_send); 1.2398 +Datum pgl_ebox_send(PG_FUNCTION_ARGS) { 1.2399 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 1.2400 + StringInfoData buf; 1.2401 + pq_begintypsend(&buf); 1.2402 + pq_sendfloat8(&buf, box->lat_min); 1.2403 + pq_sendfloat8(&buf, box->lat_max); 1.2404 + pq_sendfloat8(&buf, box->lon_min); 1.2405 + pq_sendfloat8(&buf, box->lon_max); 1.2406 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 1.2407 +} 1.2408 + 1.2409 +/* binary output function for circle ("ecircle") */ 1.2410 +PG_FUNCTION_INFO_V1(pgl_ecircle_send); 1.2411 +Datum pgl_ecircle_send(PG_FUNCTION_ARGS) { 1.2412 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 1.2413 + StringInfoData buf; 1.2414 + pq_begintypsend(&buf); 1.2415 + pq_sendfloat8(&buf, circle->center.lat); 1.2416 + pq_sendfloat8(&buf, circle->center.lon); 1.2417 + pq_sendfloat8(&buf, circle->radius); 1.2418 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 1.2419 +} 1.2420 + 1.2421 +/* TODO: binary send functions for cluster */ 1.2422 + 1.2423 +/* cast point ("epoint") to box ("ebox") */ 1.2424 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox); 1.2425 +Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) { 1.2426 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2427 + pgl_box *box = palloc(sizeof(pgl_box)); 1.2428 + box->lat_min = point->lat; 1.2429 + box->lat_max = point->lat; 1.2430 + box->lon_min = point->lon; 1.2431 + box->lon_max = point->lon; 1.2432 + PG_RETURN_POINTER(box); 1.2433 +} 1.2434 + 1.2435 +/* cast point ("epoint") to circle ("ecircle") */ 1.2436 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle); 1.2437 +Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) { 1.2438 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2439 + pgl_circle *circle = palloc(sizeof(pgl_box)); 1.2440 + circle->center = *point; 1.2441 + circle->radius = 0; 1.2442 + PG_RETURN_POINTER(circle); 1.2443 +} 1.2444 + 1.2445 +/* cast point ("epoint") to cluster ("ecluster") */ 1.2446 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster); 1.2447 +Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) { 1.2448 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2449 + pgl_newentry entry; 1.2450 + pgl_cluster *cluster; 1.2451 + entry.entrytype = PGL_ENTRY_POINT; 1.2452 + entry.npoints = 1; 1.2453 + entry.points = point; 1.2454 + cluster = pgl_new_cluster(1, &entry); 1.2455 + pgl_finalize_cluster(cluster); /* NOTE: should not fail */ 1.2456 + PG_RETURN_POINTER(cluster); 1.2457 +} 1.2458 + 1.2459 +/* cast box ("ebox") to cluster ("ecluster") */ 1.2460 +#define pgl_ebox_to_ecluster_macro(i, a, b) \ 1.2461 + entries[i].entrytype = PGL_ENTRY_POLYGON; \ 1.2462 + entries[i].npoints = 4; \ 1.2463 + entries[i].points = points[i]; \ 1.2464 + points[i][0].lat = box->lat_min; \ 1.2465 + points[i][0].lon = (a); \ 1.2466 + points[i][1].lat = box->lat_min; \ 1.2467 + points[i][1].lon = (b); \ 1.2468 + points[i][2].lat = box->lat_max; \ 1.2469 + points[i][2].lon = (b); \ 1.2470 + points[i][3].lat = box->lat_max; \ 1.2471 + points[i][3].lon = (a); 1.2472 +PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster); 1.2473 +Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) { 1.2474 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 1.2475 + double lon, dlon; 1.2476 + int nentries; 1.2477 + pgl_newentry entries[3]; 1.2478 + pgl_point points[3][4]; 1.2479 + pgl_cluster *cluster; 1.2480 + if (box->lat_min > box->lat_max) { 1.2481 + nentries = 0; 1.2482 + } else if (box->lon_min > box->lon_max) { 1.2483 + if (box->lon_min < 0) { 1.2484 + lon = pgl_round((box->lon_min + 180) / 2.0); 1.2485 + nentries = 3; 1.2486 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 1.2487 + pgl_ebox_to_ecluster_macro(1, lon, 180); 1.2488 + pgl_ebox_to_ecluster_macro(2, -180, box->lon_max); 1.2489 + } else if (box->lon_max > 0) { 1.2490 + lon = pgl_round((box->lon_max - 180) / 2.0); 1.2491 + nentries = 3; 1.2492 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 1.2493 + pgl_ebox_to_ecluster_macro(1, -180, lon); 1.2494 + pgl_ebox_to_ecluster_macro(2, lon, box->lon_max); 1.2495 + } else { 1.2496 + nentries = 2; 1.2497 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 1.2498 + pgl_ebox_to_ecluster_macro(1, -180, box->lon_max); 1.2499 + } 1.2500 + } else { 1.2501 + dlon = pgl_round(box->lon_max - box->lon_min); 1.2502 + if (dlon < 180) { 1.2503 + nentries = 1; 1.2504 + pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max); 1.2505 + } else { 1.2506 + lon = pgl_round((box->lon_min + box->lon_max) / 2.0); 1.2507 + if ( 1.2508 + pgl_round(lon - box->lon_min) < 180 && 1.2509 + pgl_round(box->lon_max - lon) < 180 1.2510 + ) { 1.2511 + nentries = 2; 1.2512 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 1.2513 + pgl_ebox_to_ecluster_macro(1, lon, box->lon_max); 1.2514 + } else { 1.2515 + nentries = 3; 1.2516 + pgl_ebox_to_ecluster_macro(0, box->lon_min, -60); 1.2517 + pgl_ebox_to_ecluster_macro(1, -60, 60); 1.2518 + pgl_ebox_to_ecluster_macro(2, 60, box->lon_max); 1.2519 + } 1.2520 + } 1.2521 + } 1.2522 + cluster = pgl_new_cluster(nentries, entries); 1.2523 + pgl_finalize_cluster(cluster); /* NOTE: should not fail */ 1.2524 + PG_RETURN_POINTER(cluster); 1.2525 +} 1.2526 + 1.2527 +/* extract latitude from point ("epoint") */ 1.2528 +PG_FUNCTION_INFO_V1(pgl_epoint_lat); 1.2529 +Datum pgl_epoint_lat(PG_FUNCTION_ARGS) { 1.2530 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat); 1.2531 +} 1.2532 + 1.2533 +/* extract longitude from point ("epoint") */ 1.2534 +PG_FUNCTION_INFO_V1(pgl_epoint_lon); 1.2535 +Datum pgl_epoint_lon(PG_FUNCTION_ARGS) { 1.2536 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon); 1.2537 +} 1.2538 + 1.2539 +/* extract minimum latitude from box ("ebox") */ 1.2540 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_min); 1.2541 +Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) { 1.2542 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min); 1.2543 +} 1.2544 + 1.2545 +/* extract maximum latitude from box ("ebox") */ 1.2546 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_max); 1.2547 +Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) { 1.2548 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max); 1.2549 +} 1.2550 + 1.2551 +/* extract minimum longitude from box ("ebox") */ 1.2552 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_min); 1.2553 +Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) { 1.2554 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min); 1.2555 +} 1.2556 + 1.2557 +/* extract maximum longitude from box ("ebox") */ 1.2558 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_max); 1.2559 +Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) { 1.2560 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max); 1.2561 +} 1.2562 + 1.2563 +/* extract center point from circle ("ecircle") */ 1.2564 +PG_FUNCTION_INFO_V1(pgl_ecircle_center); 1.2565 +Datum pgl_ecircle_center(PG_FUNCTION_ARGS) { 1.2566 + PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center)); 1.2567 +} 1.2568 + 1.2569 +/* extract radius from circle ("ecircle") */ 1.2570 +PG_FUNCTION_INFO_V1(pgl_ecircle_radius); 1.2571 +Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) { 1.2572 + PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius); 1.2573 +} 1.2574 + 1.2575 +/* check if point is inside box (overlap operator "&&") in SQL */ 1.2576 +PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap); 1.2577 +Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) { 1.2578 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2579 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1); 1.2580 + PG_RETURN_BOOL(pgl_point_in_box(point, box)); 1.2581 +} 1.2582 + 1.2583 +/* check if point is inside circle (overlap operator "&&") in SQL */ 1.2584 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap); 1.2585 +Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) { 1.2586 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2587 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 1.2588 + PG_RETURN_BOOL( 1.2589 + pgl_distance( 1.2590 + point->lat, point->lon, 1.2591 + circle->center.lat, circle->center.lon 1.2592 + ) <= circle->radius 1.2593 + ); 1.2594 +} 1.2595 + 1.2596 +/* check if point is inside cluster (overlap operator "&&") in SQL */ 1.2597 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap); 1.2598 +Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) { 1.2599 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2600 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2601 + bool retval; 1.2602 + /* points outside bounding circle are always assumed to be non-overlapping */ 1.2603 + if ( 1.2604 + pgl_distance( 1.2605 + point->lat, point->lon, 1.2606 + cluster->bounding.center.lat, cluster->bounding.center.lon 1.2607 + ) > cluster->bounding.radius 1.2608 + ) retval = false; 1.2609 + else retval = pgl_point_in_cluster(point, cluster, false); 1.2610 + PG_FREE_IF_COPY(cluster, 1); 1.2611 + PG_RETURN_BOOL(retval); 1.2612 +} 1.2613 + 1.2614 +/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */ 1.2615 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap); 1.2616 +Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) { 1.2617 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2618 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2619 + bool retval = pgl_distance( 1.2620 + point->lat, point->lon, 1.2621 + cluster->bounding.center.lat, cluster->bounding.center.lon 1.2622 + ) <= cluster->bounding.radius; 1.2623 + PG_FREE_IF_COPY(cluster, 1); 1.2624 + PG_RETURN_BOOL(retval); 1.2625 +} 1.2626 + 1.2627 +/* check if two boxes overlap (overlap operator "&&") in SQL */ 1.2628 +PG_FUNCTION_INFO_V1(pgl_ebox_overlap); 1.2629 +Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) { 1.2630 + pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0); 1.2631 + pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1); 1.2632 + PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2)); 1.2633 +} 1.2634 + 1.2635 +/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */ 1.2636 +PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap); 1.2637 +Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) { 1.2638 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 1.2639 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 1.2640 + PG_RETURN_BOOL( 1.2641 + pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius 1.2642 + ); 1.2643 +} 1.2644 + 1.2645 +/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */ 1.2646 +PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap); 1.2647 +Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) { 1.2648 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 1.2649 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2650 + bool retval = pgl_estimate_point_box_distance( 1.2651 + &cluster->bounding.center, 1.2652 + box 1.2653 + ) <= cluster->bounding.radius; 1.2654 + PG_FREE_IF_COPY(cluster, 1); 1.2655 + PG_RETURN_BOOL(retval); 1.2656 +} 1.2657 + 1.2658 +/* check if two circles overlap (overlap operator "&&") in SQL */ 1.2659 +PG_FUNCTION_INFO_V1(pgl_ecircle_overlap); 1.2660 +Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) { 1.2661 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 1.2662 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 1.2663 + PG_RETURN_BOOL( 1.2664 + pgl_distance( 1.2665 + circle1->center.lat, circle1->center.lon, 1.2666 + circle2->center.lat, circle2->center.lon 1.2667 + ) <= circle1->radius + circle2->radius 1.2668 + ); 1.2669 +} 1.2670 + 1.2671 +/* check if circle and cluster overlap (overlap operator "&&") in SQL */ 1.2672 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap); 1.2673 +Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) { 1.2674 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 1.2675 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2676 + bool retval; 1.2677 + /* points outside bounding circle (with margin for flattening) are always 1.2678 + assumed to be non-overlapping */ 1.2679 + if ( 1.2680 + (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */ 1.2681 + pgl_distance( 1.2682 + circle->center.lat, circle->center.lon, 1.2683 + cluster->bounding.center.lat, cluster->bounding.center.lon 1.2684 + ) > circle->radius + cluster->bounding.radius 1.2685 + ) retval = false; 1.2686 + else retval = pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius; 1.2687 + PG_FREE_IF_COPY(cluster, 1); 1.2688 + PG_RETURN_BOOL(retval); 1.2689 +} 1.2690 + 1.2691 +/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */ 1.2692 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap); 1.2693 +Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) { 1.2694 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 1.2695 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2696 + bool retval = ( 1.2697 + (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */ 1.2698 + pgl_distance( 1.2699 + circle->center.lat, circle->center.lon, 1.2700 + cluster->bounding.center.lat, cluster->bounding.center.lon 1.2701 + ) 1.2702 + ) <= circle->radius + cluster->bounding.radius; 1.2703 + PG_FREE_IF_COPY(cluster, 1); 1.2704 + PG_RETURN_BOOL(retval); 1.2705 +} 1.2706 + 1.2707 +/* check if two clusters overlap (overlap operator "&&") in SQL */ 1.2708 +PG_FUNCTION_INFO_V1(pgl_ecluster_overlap); 1.2709 +Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) { 1.2710 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2711 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2712 + bool retval; 1.2713 + /* clusters with non-touching bounding circles (with margin for flattening) 1.2714 + are always assumed to be non-overlapping */ 1.2715 + if ( 1.2716 + (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */ 1.2717 + pgl_distance( 1.2718 + cluster1->bounding.center.lat, cluster1->bounding.center.lon, 1.2719 + cluster2->bounding.center.lat, cluster2->bounding.center.lon 1.2720 + ) > cluster1->bounding.radius + cluster2->bounding.radius 1.2721 + ) retval = false; 1.2722 + else retval = pgl_clusters_overlap(cluster1, cluster2); 1.2723 + PG_FREE_IF_COPY(cluster1, 0); 1.2724 + PG_FREE_IF_COPY(cluster2, 1); 1.2725 + PG_RETURN_BOOL(retval); 1.2726 +} 1.2727 + 1.2728 +/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */ 1.2729 +PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap); 1.2730 +Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) { 1.2731 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2732 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2733 + bool retval = ( 1.2734 + (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */ 1.2735 + pgl_distance( 1.2736 + cluster1->bounding.center.lat, cluster1->bounding.center.lon, 1.2737 + cluster2->bounding.center.lat, cluster2->bounding.center.lon 1.2738 + ) 1.2739 + ) <= cluster1->bounding.radius + cluster2->bounding.radius; 1.2740 + PG_FREE_IF_COPY(cluster1, 0); 1.2741 + PG_FREE_IF_COPY(cluster2, 1); 1.2742 + PG_RETURN_BOOL(retval); 1.2743 +} 1.2744 + 1.2745 +/* check if second cluster is in first cluster (cont. operator "@>) in SQL */ 1.2746 +PG_FUNCTION_INFO_V1(pgl_ecluster_contains); 1.2747 +Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) { 1.2748 + pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2749 + pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2750 + bool retval; 1.2751 + /* clusters with non-touching bounding circles are always assumed to be 1.2752 + non-overlapping */ 1.2753 + if ( 1.2754 + (1.0-PGL_SPHEROID_F) * /* margin for flattening and approximation */ 1.2755 + pgl_distance( 1.2756 + outer->bounding.center.lat, outer->bounding.center.lon, 1.2757 + inner->bounding.center.lat, inner->bounding.center.lon 1.2758 + ) > outer->bounding.radius + inner->bounding.radius 1.2759 + ) retval = false; 1.2760 + else retval = pgl_cluster_in_cluster(outer, inner); 1.2761 + PG_FREE_IF_COPY(outer, 0); 1.2762 + PG_FREE_IF_COPY(inner, 1); 1.2763 + PG_RETURN_BOOL(retval); 1.2764 +} 1.2765 + 1.2766 +/* calculate distance between two points ("<->" operator) in SQL */ 1.2767 +PG_FUNCTION_INFO_V1(pgl_epoint_distance); 1.2768 +Datum pgl_epoint_distance(PG_FUNCTION_ARGS) { 1.2769 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 1.2770 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 1.2771 + PG_RETURN_FLOAT8(pgl_distance( 1.2772 + point1->lat, point1->lon, point2->lat, point2->lon 1.2773 + )); 1.2774 +} 1.2775 + 1.2776 +/* calculate point to circle distance ("<->" operator) in SQL */ 1.2777 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance); 1.2778 +Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) { 1.2779 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2780 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 1.2781 + double distance = pgl_distance( 1.2782 + point->lat, point->lon, circle->center.lat, circle->center.lon 1.2783 + ) - circle->radius; 1.2784 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 1.2785 +} 1.2786 + 1.2787 +/* calculate point to cluster distance ("<->" operator) in SQL */ 1.2788 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance); 1.2789 +Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) { 1.2790 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 1.2791 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2792 + double distance = pgl_point_cluster_distance(point, cluster); 1.2793 + PG_FREE_IF_COPY(cluster, 1); 1.2794 + PG_RETURN_FLOAT8(distance); 1.2795 +} 1.2796 + 1.2797 +/* calculate distance between two circles ("<->" operator) in SQL */ 1.2798 +PG_FUNCTION_INFO_V1(pgl_ecircle_distance); 1.2799 +Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) { 1.2800 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 1.2801 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 1.2802 + double distance = pgl_distance( 1.2803 + circle1->center.lat, circle1->center.lon, 1.2804 + circle2->center.lat, circle2->center.lon 1.2805 + ) - (circle1->radius + circle2->radius); 1.2806 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 1.2807 +} 1.2808 + 1.2809 +/* calculate circle to cluster distance ("<->" operator) in SQL */ 1.2810 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance); 1.2811 +Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) { 1.2812 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 1.2813 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2814 + double distance = ( 1.2815 + pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius 1.2816 + ); 1.2817 + PG_FREE_IF_COPY(cluster, 1); 1.2818 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 1.2819 +} 1.2820 + 1.2821 +/* calculate distance between two clusters ("<->" operator) in SQL */ 1.2822 +PG_FUNCTION_INFO_V1(pgl_ecluster_distance); 1.2823 +Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) { 1.2824 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2825 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2826 + double retval = pgl_cluster_distance(cluster1, cluster2); 1.2827 + PG_FREE_IF_COPY(cluster1, 0); 1.2828 + PG_FREE_IF_COPY(cluster2, 1); 1.2829 + PG_RETURN_FLOAT8(retval); 1.2830 +} 1.2831 + 1.2832 +/* calculate fair distance (see README) between cluster and point with 1.2833 + precision denoted by sample count ("<=> operator) in SQL */ 1.2834 +PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_sc_fair_distance); 1.2835 +Datum pgl_ecluster_epoint_sc_fair_distance(PG_FUNCTION_ARGS) { 1.2836 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1.2837 + pgl_point_sc *search = (pgl_point_sc *)PG_GETARG_POINTER(1); 1.2838 + double retval = pgl_fair_distance( 1.2839 + &search->point, cluster, search->samples 1.2840 + ); 1.2841 + PG_FREE_IF_COPY(cluster, 0); 1.2842 + PG_RETURN_FLOAT8(retval); 1.2843 +} 1.2844 + 1.2845 + 1.2846 +/*-----------------------------------------------------------* 1.2847 + * B-tree comparison operators and index support functions * 1.2848 + *-----------------------------------------------------------*/ 1.2849 + 1.2850 +/* macro for a B-tree operator (without detoasting) */ 1.2851 +#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \ 1.2852 + PG_FUNCTION_INFO_V1(func); \ 1.2853 + Datum func(PG_FUNCTION_ARGS) { \ 1.2854 + type *a = (type *)PG_GETARG_POINTER(0); \ 1.2855 + type *b = (type *)PG_GETARG_POINTER(1); \ 1.2856 + PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \ 1.2857 + } 1.2858 + 1.2859 +/* macro for a B-tree comparison function (without detoasting) */ 1.2860 +#define PGL_BTREE_CMP(func, type, cmpfunc) \ 1.2861 + PG_FUNCTION_INFO_V1(func); \ 1.2862 + Datum func(PG_FUNCTION_ARGS) { \ 1.2863 + type *a = (type *)PG_GETARG_POINTER(0); \ 1.2864 + type *b = (type *)PG_GETARG_POINTER(1); \ 1.2865 + PG_RETURN_INT32(cmpfunc(a, b)); \ 1.2866 + } 1.2867 + 1.2868 +/* macro for a B-tree operator (with detoasting) */ 1.2869 +#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \ 1.2870 + PG_FUNCTION_INFO_V1(func); \ 1.2871 + Datum func(PG_FUNCTION_ARGS) { \ 1.2872 + bool res; \ 1.2873 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 1.2874 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 1.2875 + res = cmpfunc(a, b) oper 0; \ 1.2876 + PG_FREE_IF_COPY(a, 0); \ 1.2877 + PG_FREE_IF_COPY(b, 1); \ 1.2878 + PG_RETURN_BOOL(res); \ 1.2879 + } 1.2880 + 1.2881 +/* macro for a B-tree comparison function (with detoasting) */ 1.2882 +#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \ 1.2883 + PG_FUNCTION_INFO_V1(func); \ 1.2884 + Datum func(PG_FUNCTION_ARGS) { \ 1.2885 + int32_t res; \ 1.2886 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 1.2887 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 1.2888 + res = cmpfunc(a, b); \ 1.2889 + PG_FREE_IF_COPY(a, 0); \ 1.2890 + PG_FREE_IF_COPY(b, 1); \ 1.2891 + PG_RETURN_INT32(res); \ 1.2892 + } 1.2893 + 1.2894 +/* B-tree operators and comparison function for point */ 1.2895 +PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <) 1.2896 +PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=) 1.2897 +PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==) 1.2898 +PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=) 1.2899 +PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=) 1.2900 +PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >) 1.2901 +PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp) 1.2902 + 1.2903 +/* B-tree operators and comparison function for box */ 1.2904 +PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <) 1.2905 +PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=) 1.2906 +PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==) 1.2907 +PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=) 1.2908 +PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=) 1.2909 +PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >) 1.2910 +PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp) 1.2911 + 1.2912 +/* B-tree operators and comparison function for circle */ 1.2913 +PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <) 1.2914 +PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=) 1.2915 +PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==) 1.2916 +PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=) 1.2917 +PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=) 1.2918 +PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >) 1.2919 +PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp) 1.2920 + 1.2921 + 1.2922 +/*--------------------------------* 1.2923 + * GiST index support functions * 1.2924 + *--------------------------------*/ 1.2925 + 1.2926 +/* GiST "consistent" support function */ 1.2927 +PG_FUNCTION_INFO_V1(pgl_gist_consistent); 1.2928 +Datum pgl_gist_consistent(PG_FUNCTION_ARGS) { 1.2929 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 1.2930 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 1.2931 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 1.2932 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 1.2933 + /* demand recheck because index and query methods are lossy */ 1.2934 + *recheck = true; 1.2935 + /* strategy number aliases for different operators using the same strategy */ 1.2936 + strategy %= 100; 1.2937 + /* strategy number 11: equality of two points */ 1.2938 + if (strategy == 11) { 1.2939 + /* query datum is another point */ 1.2940 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 1.2941 + /* convert other point to key */ 1.2942 + pgl_pointkey querykey; 1.2943 + pgl_point_to_key(query, querykey); 1.2944 + /* return true if both keys overlap */ 1.2945 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 1.2946 + } 1.2947 + /* strategy number 13: equality of two circles */ 1.2948 + if (strategy == 13) { 1.2949 + /* query datum is another circle */ 1.2950 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 1.2951 + /* convert other circle to key */ 1.2952 + pgl_areakey querykey; 1.2953 + pgl_circle_to_key(query, querykey); 1.2954 + /* return true if both keys overlap */ 1.2955 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 1.2956 + } 1.2957 + /* for all remaining strategies, keys on empty objects produce no match */ 1.2958 + /* (check necessary because query radius may be infinite) */ 1.2959 + if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false); 1.2960 + /* strategy number 21: overlapping with point */ 1.2961 + if (strategy == 21) { 1.2962 + /* query datum is a point */ 1.2963 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 1.2964 + /* return true if estimated distance (allowed to be smaller than real 1.2965 + distance) between index key and point is zero */ 1.2966 + PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0); 1.2967 + } 1.2968 + /* strategy number 22: (point) overlapping with box */ 1.2969 + if (strategy == 22) { 1.2970 + /* query datum is a box */ 1.2971 + pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1); 1.2972 + /* determine bounding box of indexed key */ 1.2973 + pgl_box keybox; 1.2974 + pgl_key_to_box(key, &keybox); 1.2975 + /* return true if query box overlaps with bounding box of indexed key */ 1.2976 + PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox)); 1.2977 + } 1.2978 + /* strategy number 23: overlapping with circle */ 1.2979 + if (strategy == 23) { 1.2980 + /* query datum is a circle */ 1.2981 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 1.2982 + /* return true if estimated distance (allowed to be smaller than real 1.2983 + distance) between index key and circle center is smaller than radius */ 1.2984 + PG_RETURN_BOOL( 1.2985 + (1.0-PGL_SPHEROID_F) * /* safety margin for lossy operator */ 1.2986 + pgl_estimate_key_distance(key, &(query->center)) 1.2987 + <= query->radius 1.2988 + ); 1.2989 + } 1.2990 + /* strategy number 24: overlapping with cluster */ 1.2991 + if (strategy == 24) { 1.2992 + bool retval; /* return value */ 1.2993 + /* query datum is a cluster */ 1.2994 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.2995 + /* return true if estimated distance (allowed to be smaller than real 1.2996 + distance) between index key and circle center is smaller than radius */ 1.2997 + retval = ( 1.2998 + (1.0-PGL_SPHEROID_F) * /* safety margin for lossy operator */ 1.2999 + pgl_estimate_key_distance(key, &(query->bounding.center)) 1.3000 + <= query->bounding.radius 1.3001 + ); 1.3002 + PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */ 1.3003 + PG_RETURN_BOOL(retval); 1.3004 + } 1.3005 + /* throw error for any unknown strategy number */ 1.3006 + elog(ERROR, "unrecognized strategy number: %d", strategy); 1.3007 +} 1.3008 + 1.3009 +/* GiST "union" support function */ 1.3010 +PG_FUNCTION_INFO_V1(pgl_gist_union); 1.3011 +Datum pgl_gist_union(PG_FUNCTION_ARGS) { 1.3012 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 1.3013 + pgl_keyptr out; /* return value (to be palloc'ed) */ 1.3014 + int i; 1.3015 + /* determine key size */ 1.3016 + size_t keysize = PGL_KEY_IS_AREAKEY( 1.3017 + (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key) 1.3018 + ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey); 1.3019 + /* begin with first key as result */ 1.3020 + out = palloc(keysize); 1.3021 + memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize); 1.3022 + /* unite current result with second, third, etc. key */ 1.3023 + for (i=1; i<entryvec->n; i++) { 1.3024 + pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)); 1.3025 + } 1.3026 + /* return result */ 1.3027 + PG_RETURN_POINTER(out); 1.3028 +} 1.3029 + 1.3030 +/* GiST "compress" support function for indicis on points */ 1.3031 +PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint); 1.3032 +Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) { 1.3033 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 1.3034 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 1.3035 + /* only transform new leaves */ 1.3036 + if (entry->leafkey) { 1.3037 + /* get point to be transformed */ 1.3038 + pgl_point *point = (pgl_point *)DatumGetPointer(entry->key); 1.3039 + /* allocate memory for key */ 1.3040 + pgl_keyptr key = palloc(sizeof(pgl_pointkey)); 1.3041 + /* transform point to key */ 1.3042 + pgl_point_to_key(point, key); 1.3043 + /* create new GISTENTRY structure as return value */ 1.3044 + retval = palloc(sizeof(GISTENTRY)); 1.3045 + gistentryinit( 1.3046 + *retval, PointerGetDatum(key), 1.3047 + entry->rel, entry->page, entry->offset, false 1.3048 + ); 1.3049 + } else { 1.3050 + /* inner nodes have already been transformed */ 1.3051 + retval = entry; 1.3052 + } 1.3053 + /* return pointer to old or new GISTENTRY structure */ 1.3054 + PG_RETURN_POINTER(retval); 1.3055 +} 1.3056 + 1.3057 +/* GiST "compress" support function for indicis on circles */ 1.3058 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle); 1.3059 +Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) { 1.3060 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 1.3061 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 1.3062 + /* only transform new leaves */ 1.3063 + if (entry->leafkey) { 1.3064 + /* get circle to be transformed */ 1.3065 + pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key); 1.3066 + /* allocate memory for key */ 1.3067 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 1.3068 + /* transform circle to key */ 1.3069 + pgl_circle_to_key(circle, key); 1.3070 + /* create new GISTENTRY structure as return value */ 1.3071 + retval = palloc(sizeof(GISTENTRY)); 1.3072 + gistentryinit( 1.3073 + *retval, PointerGetDatum(key), 1.3074 + entry->rel, entry->page, entry->offset, false 1.3075 + ); 1.3076 + } else { 1.3077 + /* inner nodes have already been transformed */ 1.3078 + retval = entry; 1.3079 + } 1.3080 + /* return pointer to old or new GISTENTRY structure */ 1.3081 + PG_RETURN_POINTER(retval); 1.3082 +} 1.3083 + 1.3084 +/* GiST "compress" support function for indices on clusters */ 1.3085 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster); 1.3086 +Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) { 1.3087 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 1.3088 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 1.3089 + /* only transform new leaves */ 1.3090 + if (entry->leafkey) { 1.3091 + /* get cluster to be transformed (detoasting necessary!) */ 1.3092 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key); 1.3093 + /* allocate memory for key */ 1.3094 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 1.3095 + /* transform cluster to key */ 1.3096 + pgl_circle_to_key(&(cluster->bounding), key); 1.3097 + /* create new GISTENTRY structure as return value */ 1.3098 + retval = palloc(sizeof(GISTENTRY)); 1.3099 + gistentryinit( 1.3100 + *retval, PointerGetDatum(key), 1.3101 + entry->rel, entry->page, entry->offset, false 1.3102 + ); 1.3103 + /* free detoasted datum */ 1.3104 + if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster); 1.3105 + } else { 1.3106 + /* inner nodes have already been transformed */ 1.3107 + retval = entry; 1.3108 + } 1.3109 + /* return pointer to old or new GISTENTRY structure */ 1.3110 + PG_RETURN_POINTER(retval); 1.3111 +} 1.3112 + 1.3113 +/* GiST "decompress" support function for indices */ 1.3114 +PG_FUNCTION_INFO_V1(pgl_gist_decompress); 1.3115 +Datum pgl_gist_decompress(PG_FUNCTION_ARGS) { 1.3116 + /* return passed pointer without transformation */ 1.3117 + PG_RETURN_POINTER(PG_GETARG_POINTER(0)); 1.3118 +} 1.3119 + 1.3120 +/* GiST "penalty" support function */ 1.3121 +PG_FUNCTION_INFO_V1(pgl_gist_penalty); 1.3122 +Datum pgl_gist_penalty(PG_FUNCTION_ARGS) { 1.3123 + GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0); 1.3124 + GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1); 1.3125 + float *penalty = (float *)PG_GETARG_POINTER(2); 1.3126 + /* get original key and key to insert */ 1.3127 + pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key); 1.3128 + pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key); 1.3129 + /* copy original key */ 1.3130 + union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key; 1.3131 + if (PGL_KEY_IS_AREAKEY(orig)) { 1.3132 + memcpy(union_key.areakey, orig, sizeof(union_key.areakey)); 1.3133 + } else { 1.3134 + memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey)); 1.3135 + } 1.3136 + /* calculate union of both keys */ 1.3137 + pgl_unite_keys((pgl_keyptr)&union_key, new); 1.3138 + /* penalty equal to reduction of key length (logarithm of added area) */ 1.3139 + /* (return value by setting referenced value and returning pointer) */ 1.3140 + *penalty = ( 1.3141 + PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key) 1.3142 + ); 1.3143 + PG_RETURN_POINTER(penalty); 1.3144 +} 1.3145 + 1.3146 +/* GiST "picksplit" support function */ 1.3147 +PG_FUNCTION_INFO_V1(pgl_gist_picksplit); 1.3148 +Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) { 1.3149 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 1.3150 + GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1); 1.3151 + OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (exclusive) */ 1.3152 + union { 1.3153 + pgl_pointkey pointkey; 1.3154 + pgl_areakey areakey; 1.3155 + } union_all; /* union of all keys (to be calculated from scratch) 1.3156 + (later cut in half) */ 1.3157 + int is_areakey = PGL_KEY_IS_AREAKEY( 1.3158 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key) 1.3159 + ); 1.3160 + int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey); 1.3161 + pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */ 1.3162 + pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */ 1.3163 + pgl_keyptr key; /* current key to be processed */ 1.3164 + /* allocate memory for array of left and right keys, set counts to zero */ 1.3165 + v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 1.3166 + v->spl_nleft = 0; 1.3167 + v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 1.3168 + v->spl_nright = 0; 1.3169 + /* calculate union of all keys from scratch */ 1.3170 + memcpy( 1.3171 + (pgl_keyptr)&union_all, 1.3172 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key), 1.3173 + keysize 1.3174 + ); 1.3175 + for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) { 1.3176 + pgl_unite_keys( 1.3177 + (pgl_keyptr)&union_all, 1.3178 + (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key) 1.3179 + ); 1.3180 + } 1.3181 + /* check if trivial split is necessary due to exhausted key length */ 1.3182 + /* (Note: keys for empty objects must have node depth set to maximum) */ 1.3183 + if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == ( 1.3184 + is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH 1.3185 + )) { 1.3186 + /* half of all keys go left */ 1.3187 + for ( 1.3188 + i=FirstOffsetNumber; 1.3189 + i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 1.3190 + i=OffsetNumberNext(i) 1.3191 + ) { 1.3192 + /* pointer to current key */ 1.3193 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 1.3194 + /* update unionL */ 1.3195 + /* check if key is first key that goes left */ 1.3196 + if (!v->spl_nleft) { 1.3197 + /* first key that goes left is just copied to unionL */ 1.3198 + memcpy(unionL, key, keysize); 1.3199 + } else { 1.3200 + /* unite current value and next key */ 1.3201 + pgl_unite_keys(unionL, key); 1.3202 + } 1.3203 + /* append offset number to list of keys that go left */ 1.3204 + v->spl_left[v->spl_nleft++] = i; 1.3205 + } 1.3206 + /* other half goes right */ 1.3207 + for ( 1.3208 + i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 1.3209 + i<entryvec->n; 1.3210 + i=OffsetNumberNext(i) 1.3211 + ) { 1.3212 + /* pointer to current key */ 1.3213 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 1.3214 + /* update unionR */ 1.3215 + /* check if key is first key that goes right */ 1.3216 + if (!v->spl_nright) { 1.3217 + /* first key that goes right is just copied to unionR */ 1.3218 + memcpy(unionR, key, keysize); 1.3219 + } else { 1.3220 + /* unite current value and next key */ 1.3221 + pgl_unite_keys(unionR, key); 1.3222 + } 1.3223 + /* append offset number to list of keys that go right */ 1.3224 + v->spl_right[v->spl_nright++] = i; 1.3225 + } 1.3226 + } 1.3227 + /* otherwise, a non-trivial split is possible */ 1.3228 + else { 1.3229 + /* cut covered area in half */ 1.3230 + /* (union_all then refers to area of keys that go left) */ 1.3231 + /* check if union of all keys covers empty and non-empty objects */ 1.3232 + if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) { 1.3233 + /* if yes, split into empty and non-empty objects */ 1.3234 + pgl_key_set_empty((pgl_keyptr)&union_all); 1.3235 + } else { 1.3236 + /* otherwise split by next bit */ 1.3237 + ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++; 1.3238 + /* NOTE: type bit conserved */ 1.3239 + } 1.3240 + /* determine for each key if it goes left or right */ 1.3241 + for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) { 1.3242 + /* pointer to current key */ 1.3243 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 1.3244 + /* keys within one half of the area go left */ 1.3245 + if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) { 1.3246 + /* update unionL */ 1.3247 + /* check if key is first key that goes left */ 1.3248 + if (!v->spl_nleft) { 1.3249 + /* first key that goes left is just copied to unionL */ 1.3250 + memcpy(unionL, key, keysize); 1.3251 + } else { 1.3252 + /* unite current value of unionL and processed key */ 1.3253 + pgl_unite_keys(unionL, key); 1.3254 + } 1.3255 + /* append offset number to list of keys that go left */ 1.3256 + v->spl_left[v->spl_nleft++] = i; 1.3257 + } 1.3258 + /* the other keys go right */ 1.3259 + else { 1.3260 + /* update unionR */ 1.3261 + /* check if key is first key that goes right */ 1.3262 + if (!v->spl_nright) { 1.3263 + /* first key that goes right is just copied to unionR */ 1.3264 + memcpy(unionR, key, keysize); 1.3265 + } else { 1.3266 + /* unite current value of unionR and processed key */ 1.3267 + pgl_unite_keys(unionR, key); 1.3268 + } 1.3269 + /* append offset number to list of keys that go right */ 1.3270 + v->spl_right[v->spl_nright++] = i; 1.3271 + } 1.3272 + } 1.3273 + } 1.3274 + /* store unions in return value */ 1.3275 + v->spl_ldatum = PointerGetDatum(unionL); 1.3276 + v->spl_rdatum = PointerGetDatum(unionR); 1.3277 + /* return all results */ 1.3278 + PG_RETURN_POINTER(v); 1.3279 +} 1.3280 + 1.3281 +/* GiST "same"/"equal" support function */ 1.3282 +PG_FUNCTION_INFO_V1(pgl_gist_same); 1.3283 +Datum pgl_gist_same(PG_FUNCTION_ARGS) { 1.3284 + pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0); 1.3285 + pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1); 1.3286 + bool *boolptr = (bool *)PG_GETARG_POINTER(2); 1.3287 + /* two keys are equal if they are binary equal */ 1.3288 + /* (return result by setting referenced boolean and returning pointer) */ 1.3289 + *boolptr = !memcmp( 1.3290 + key1, 1.3291 + key2, 1.3292 + PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey) 1.3293 + ); 1.3294 + PG_RETURN_POINTER(boolptr); 1.3295 +} 1.3296 + 1.3297 +/* GiST "distance" support function */ 1.3298 +PG_FUNCTION_INFO_V1(pgl_gist_distance); 1.3299 +Datum pgl_gist_distance(PG_FUNCTION_ARGS) { 1.3300 + GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0); 1.3301 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 1.3302 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 1.3303 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 1.3304 + double distance; /* return value */ 1.3305 + /* demand recheck because distance is just an estimation */ 1.3306 + /* (real distance may be bigger) */ 1.3307 + *recheck = true; 1.3308 + /* strategy number aliases for different operators using the same strategy */ 1.3309 + strategy %= 100; 1.3310 + /* strategy number 31: distance to point */ 1.3311 + if (strategy == 31) { 1.3312 + /* query datum is a point */ 1.3313 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 1.3314 + /* use pgl_estimate_pointkey_distance() function to compute result */ 1.3315 + distance = pgl_estimate_key_distance(key, query); 1.3316 + /* avoid infinity (reserved!) */ 1.3317 + if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 1.3318 + /* return result */ 1.3319 + PG_RETURN_FLOAT8(distance); 1.3320 + } 1.3321 + /* strategy number 33: distance to circle */ 1.3322 + if (strategy == 33) { 1.3323 + /* query datum is a circle */ 1.3324 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 1.3325 + /* estimate distance to circle center and substract circle radius */ 1.3326 + distance = ( 1.3327 + pgl_estimate_key_distance(key, &(query->center)) - query->radius 1.3328 + ); 1.3329 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 1.3330 + if (distance <= 0) distance = 0; 1.3331 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 1.3332 + /* return result */ 1.3333 + PG_RETURN_FLOAT8(distance); 1.3334 + } 1.3335 + /* strategy number 34: distance to cluster */ 1.3336 + if (strategy == 34) { 1.3337 + /* query datum is a cluster */ 1.3338 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1.3339 + /* estimate distance to bounding center and substract bounding radius */ 1.3340 + distance = ( 1.3341 + pgl_estimate_key_distance(key, &(query->bounding.center)) - 1.3342 + query->bounding.radius 1.3343 + ); 1.3344 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 1.3345 + if (distance <= 0) distance = 0; 1.3346 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 1.3347 + /* free detoasted cluster (if copy) */ 1.3348 + PG_FREE_IF_COPY(query, 1); 1.3349 + /* return result */ 1.3350 + PG_RETURN_FLOAT8(distance); 1.3351 + } 1.3352 + /* throw error for any unknown strategy number */ 1.3353 + elog(ERROR, "unrecognized strategy number: %d", strategy); 1.3354 +} 1.3355 +