pgLatLon: latlon-v0010.c annotate

pgLatLon

annotate latlon-v0010.c @ 79:16787a19a325

Do not return CSTRING to static memory

Avoids wrong pfree under certain circumstances
(e.g. array_agg on empty EBOX or ECLUSTER).

author	jbe
date	Thu Oct 23 10:16:01 2025 +0200 (3 days ago)
parents	8fa0ef5e4ee6
children	b9cdc74a90db

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