jbe@0: jbe@0: /*-------------* jbe@0: * C prelude * jbe@0: *-------------*/ jbe@0: jbe@0: #include "postgres.h" jbe@0: #include "fmgr.h" jbe@0: #include "libpq/pqformat.h" jbe@0: #include "access/gist.h" jbe@0: #include "access/stratnum.h" jbe@0: #include "utils/array.h" jbe@0: #include jbe@0: jbe@0: #ifdef PG_MODULE_MAGIC jbe@0: PG_MODULE_MAGIC; jbe@0: #endif jbe@0: jbe@0: #if INT_MAX < 2147483647 jbe@0: #error Expected int type to be at least 32 bit wide jbe@0: #endif jbe@0: jbe@0: jbe@0: /*---------------------------------* jbe@0: * distance calculation on earth * jbe@0: * (using WGS-84 spheroid) * jbe@0: *---------------------------------*/ jbe@0: jbe@0: /* WGS-84 spheroid with following parameters: jbe@0: semi-major axis a = 6378137 jbe@0: semi-minor axis b = a * (1 - 1/298.257223563) jbe@0: estimated diameter = 2 * (2*a+b)/3 jbe@0: */ jbe@0: #define PGL_SPHEROID_A 6378137.0 /* semi major axis */ jbe@0: #define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */ jbe@0: #define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F)) jbe@0: #define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \ jbe@0: PGL_SPHEROID_B * PGL_SPHEROID_B ) / \ jbe@0: ( PGL_SPHEROID_A * PGL_SPHEROID_A ) ) jbe@0: #define PGL_SUBEPS2 (1.0-PGL_EPS2) jbe@0: #define PGL_DIAMETER ((4.0*PGL_SPHEROID_A + 2.0*PGL_SPHEROID_B) / 3.0) jbe@0: #define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */ jbe@0: #define PGL_FADELIMIT (PGL_DIAMETER * M_PI / 6.0) /* 1/6 circumference */ jbe@0: #define PGL_MAXDIST (PGL_DIAMETER * M_PI / 2.0) /* maximum distance */ jbe@0: jbe@0: /* calculate distance between two points on earth (given in degrees) */ jbe@0: static inline double pgl_distance( jbe@0: double lat1, double lon1, double lat2, double lon2 jbe@0: ) { jbe@0: float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin; jbe@0: float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t; jbe@0: /* normalize delta longitude (lon2 > 0 && lon1 = 0) */ jbe@0: /* lon1 = 0 (not used anymore) */ jbe@0: lon2 = fabs(lon2-lon1); jbe@0: /* convert to radians (first divide, then multiply) */ jbe@0: lat1 = (lat1 / 180.0) * M_PI; jbe@0: lat2 = (lat2 / 180.0) * M_PI; jbe@0: lon2 = (lon2 / 180.0) * M_PI; jbe@0: /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point jbe@0: operations (lon2 >= lon1 is already ensured in a previous step) */ jbe@0: if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; } jbe@0: /* calculate 3d coordinates on scaled ellipsoid which has an average diameter jbe@0: of 1.0 */ jbe@0: lat1cos = cos(lat1); lat1sin = sin(lat1); jbe@0: lat2cos = cos(lat2); lat2sin = sin(lat2); jbe@0: lon2cos = cos(lon2); lon2sin = sin(lon2); jbe@0: nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin); jbe@0: nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin); jbe@0: x1 = nphi1 * lat1cos; jbe@0: z1 = nphi1 * PGL_SUBEPS2 * lat1sin; jbe@0: x2 = nphi2 * lat2cos * lon2cos; jbe@0: y2 = nphi2 * lat2cos * lon2sin; jbe@0: z2 = nphi2 * PGL_SUBEPS2 * lat2sin; jbe@0: /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */ jbe@0: g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1)); jbe@0: /* convert tunnel distance through scaled ellipsoid to approximated surface jbe@0: distance on original ellipsoid */ jbe@0: if (g > 1.0) g = 1.0; jbe@0: s = PGL_DIAMETER * asin(g); jbe@0: /* return result only if small enough to be precise (less than 1/3 of jbe@0: maximum possible distance) */ jbe@0: if (s <= PGL_FADELIMIT) return s; jbe@0: /* calculate tunnel distance to antipodal point through scaled ellipsoid */ jbe@3: g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1)); jbe@0: /* convert tunnel distance to antipodal point through scaled ellipsoid to jbe@0: approximated surface distance to antipodal point on original ellipsoid */ jbe@0: if (g > 1.0) g = 1.0; jbe@0: t = PGL_DIAMETER * asin(g); jbe@0: /* surface distance between original points can now be approximated by jbe@0: substracting antipodal distance from maximum possible distance; jbe@0: return result only if small enough (less than 1/3 of maximum possible jbe@0: distance) */ jbe@0: if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t; jbe@0: /* otherwise crossfade direct and antipodal result to ensure monotonicity */ jbe@0: return ( jbe@0: (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) / jbe@0: (s + t - 2*PGL_FADELIMIT) jbe@0: ); jbe@0: } jbe@0: jbe@0: /* finite distance that can not be reached on earth */ jbe@0: #define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST) jbe@0: jbe@0: jbe@0: /*--------------------------------* jbe@0: * simple geographic data types * jbe@0: *--------------------------------*/ jbe@0: jbe@0: /* point on earth given by latitude and longitude in degrees */ jbe@0: /* (type "epoint" in SQL) */ jbe@0: typedef struct { jbe@0: double lat; /* between -90 and 90 (both inclusive) */ jbe@0: double lon; /* between -180 and 180 (both inclusive) */ jbe@0: } pgl_point; jbe@0: jbe@0: /* box delimited by two parallels and two meridians (all in degrees) */ jbe@0: /* (type "ebox" in SQL) */ jbe@0: typedef struct { jbe@0: double lat_min; /* between -90 and 90 (both inclusive) */ jbe@0: double lat_max; /* between -90 and 90 (both inclusive) */ jbe@0: double lon_min; /* between -180 and 180 (both inclusive) */ jbe@0: double lon_max; /* between -180 and 180 (both inclusive) */ jbe@0: /* if lat_min > lat_max, then box is empty */ jbe@0: /* if lon_min > lon_max, then 180th meridian is crossed */ jbe@0: } pgl_box; jbe@0: jbe@0: /* circle on earth surface (for radial searches with fixed radius) */ jbe@0: /* (type "ecircle" in SQL) */ jbe@0: typedef struct { jbe@0: pgl_point center; jbe@0: double radius; /* positive (including +0 but excluding -0), or -INFINITY */ jbe@0: /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point), jbe@0: zero radius (0) denotes a single point, jbe@0: a finite radius (0 < radius < INFINITY) denotes a filled circle, and jbe@0: a radius of INFINITY is valid and means complete coverage of earth. */ jbe@0: } pgl_circle; jbe@0: jbe@0: jbe@0: /*----------------------------------* jbe@0: * geographic "cluster" data type * jbe@0: *----------------------------------*/ jbe@0: jbe@0: /* A cluster is a collection of points, paths, outlines, and polygons. If two jbe@0: polygons in a cluster overlap, the area covered by both polygons does not jbe@0: belong to the cluster. This way, a cluster can be used to describe complex jbe@0: shapes like polygons with holes. Outlines are non-filled polygons. Paths are jbe@0: open by default (i.e. the last point in the list is not connected with the jbe@0: first point in the list). Note that each outline or polygon in a cluster jbe@0: must cover a longitude range of less than 180 degrees to avoid ambiguities. jbe@0: Areas which are larger may be split into multiple polygons. */ jbe@0: jbe@0: /* maximum number of points in a cluster */ jbe@0: /* (limited to avoid integer overflows, e.g. when allocating memory) */ jbe@0: #define PGL_CLUSTER_MAXPOINTS 16777216 jbe@0: jbe@0: /* types of cluster entries */ jbe@0: #define PGL_ENTRY_POINT 1 /* a point */ jbe@0: #define PGL_ENTRY_PATH 2 /* a path from first point to last point */ jbe@0: #define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */ jbe@0: #define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */ jbe@0: jbe@0: /* Entries of a cluster are described by two different structs: pgl_newentry jbe@0: and pgl_entry. The first is used only during construction of a cluster, the jbe@0: second is used in all other cases (e.g. when reading clusters from the jbe@0: database, performing operations, etc). */ jbe@0: jbe@0: /* entry for new geographic cluster during construction of that cluster */ jbe@0: typedef struct { jbe@0: int32_t entrytype; jbe@0: int32_t npoints; jbe@0: pgl_point *points; /* pointer to an array of points (pgl_point) */ jbe@0: } pgl_newentry; jbe@0: jbe@0: /* entry of geographic cluster */ jbe@0: typedef struct { jbe@0: int32_t entrytype; /* type of entry: point, path, outline, polygon */ jbe@0: int32_t npoints; /* number of stored points (set to 1 for point entry) */ jbe@0: int32_t offset; /* offset of pgl_point array from cluster base address */ jbe@0: /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */ jbe@0: } pgl_entry; jbe@0: jbe@0: /* geographic cluster which is a collection of points, (open) paths, polygons, jbe@0: and outlines (non-filled polygons) */ jbe@0: typedef struct { jbe@0: char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */ jbe@0: int32_t nentries; /* number of stored points */ jbe@0: pgl_circle bounding; /* bounding circle */ jbe@0: /* Note: bounding circle ensures alignment of pgl_cluster for points */ jbe@0: pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */ jbe@0: } pgl_cluster; jbe@0: jbe@0: /* macro to determine memory alignment of points */ jbe@0: /* (needed to store pgl_point array after entries in pgl_cluster) */ jbe@0: typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment; jbe@0: #define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned) jbe@0: jbe@0: /* macro to extract a pointer to the array of points of a cluster entry */ jbe@0: #define PGL_ENTRY_POINTS(cluster, idx) \ jbe@0: ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset)) jbe@0: jbe@0: /* convert pgl_newentry array to pgl_cluster */ jbe@0: static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) { jbe@0: int i; /* index of current entry */ jbe@0: int npoints = 0; /* number of points in whole cluster */ jbe@0: int entry_npoints; /* number of points in current entry */ jbe@0: int points_offset = PGL_POINT_ALIGNMENT * ( jbe@0: ( offsetof(pgl_cluster, entries) + jbe@0: nentries * sizeof(pgl_entry) + jbe@0: PGL_POINT_ALIGNMENT - 1 jbe@0: ) / PGL_POINT_ALIGNMENT jbe@0: ); /* offset of pgl_point array from base address (considering alignment) */ jbe@0: pgl_cluster *cluster; /* new cluster to be returned */ jbe@0: /* determine total number of points */ jbe@0: for (i=0; ientries[i].entrytype = entries[i].entrytype; jbe@0: cluster->entries[i].npoints = entry_npoints; jbe@0: /* calculate offset (in bytes) of pgl_point array */ jbe@0: cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point); jbe@0: /* copy points */ jbe@0: memcpy( jbe@0: PGL_ENTRY_POINTS(cluster, i), jbe@0: entries[i].points, jbe@0: entry_npoints * sizeof(pgl_point) jbe@0: ); jbe@0: /* update total number of points processed */ jbe@0: npoints += entry_npoints; jbe@0: } jbe@0: /* set number of entries in cluster */ jbe@0: cluster->nentries = nentries; jbe@0: /* set PostgreSQL header for variable sized data */ jbe@0: SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point)); jbe@0: /* return newly created cluster */ jbe@0: return cluster; jbe@0: } jbe@0: jbe@0: jbe@0: /*----------------------------------------* jbe@0: * C functions on geographic data types * jbe@0: *----------------------------------------*/ jbe@0: jbe@0: /* round latitude or longitude to 12 digits after decimal point */ jbe@0: static inline double pgl_round(double val) { jbe@0: return round(val * 1e12) / 1e12; jbe@0: } jbe@0: jbe@0: /* compare two points */ jbe@0: /* (equality when same point on earth is described, otherwise an arbitrary jbe@0: linear order) */ jbe@0: static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) { jbe@0: double lon1, lon2; /* modified longitudes for special cases */ jbe@0: /* use latitude as first ordering criterion */ jbe@0: if (point1->lat < point2->lat) return -1; jbe@0: if (point1->lat > point2->lat) return 1; jbe@0: /* determine modified longitudes (considering special case of poles and jbe@0: 180th meridian which can be described as W180 or E180) */ jbe@0: if (point1->lat == -90 || point1->lat == 90) lon1 = 0; jbe@0: else if (point1->lon == 180) lon1 = -180; jbe@0: else lon1 = point1->lon; jbe@0: if (point2->lat == -90 || point2->lat == 90) lon2 = 0; jbe@0: else if (point2->lon == 180) lon2 = -180; jbe@0: else lon2 = point2->lon; jbe@0: /* use (modified) longitude as secondary ordering criterion */ jbe@0: if (lon1 < lon2) return -1; jbe@0: if (lon1 > lon2) return 1; jbe@0: /* no difference found, points are equal */ jbe@0: return 0; jbe@0: } jbe@0: jbe@0: /* compare two boxes */ jbe@0: /* (equality when same box on earth is described, otherwise an arbitrary linear jbe@0: order) */ jbe@0: static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) { jbe@0: /* two empty boxes are equal, and an empty box is always considered "less jbe@0: than" a non-empty box */ jbe@0: if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1; jbe@0: if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0; jbe@0: if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1; jbe@0: /* use southern border as first ordering criterion */ jbe@0: if (box1->lat_min < box2->lat_min) return -1; jbe@0: if (box1->lat_min > box2->lat_min) return 1; jbe@0: /* use northern border as second ordering criterion */ jbe@0: if (box1->lat_max < box2->lat_max) return -1; jbe@0: if (box1->lat_max > box2->lat_max) return 1; jbe@0: /* use western border as third ordering criterion */ jbe@0: if (box1->lon_min < box2->lon_min) return -1; jbe@0: if (box1->lon_min > box2->lon_min) return 1; jbe@0: /* use eastern border as fourth ordering criterion */ jbe@0: if (box1->lon_max < box2->lon_max) return -1; jbe@0: if (box1->lon_max > box2->lon_max) return 1; jbe@0: /* no difference found, boxes are equal */ jbe@0: return 0; jbe@0: } jbe@0: jbe@0: /* compare two circles */ jbe@0: /* (equality when same circle on earth is described, otherwise an arbitrary jbe@0: linear order) */ jbe@0: static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) { jbe@0: /* two circles with same infinite radius (positive or negative infinity) are jbe@0: considered equal independently of center point */ jbe@0: if ( jbe@0: !isfinite(circle1->radius) && !isfinite(circle2->radius) && jbe@0: circle1->radius == circle2->radius jbe@0: ) return 0; jbe@0: /* use radius as first ordering criterion */ jbe@0: if (circle1->radius < circle2->radius) return -1; jbe@0: if (circle1->radius > circle2->radius) return 1; jbe@0: /* use center point as secondary ordering criterion */ jbe@0: return pgl_point_cmp(&(circle1->center), &(circle2->center)); jbe@0: } jbe@0: jbe@0: /* set box to empty box*/ jbe@0: static void pgl_box_set_empty(pgl_box *box) { jbe@0: box->lat_min = INFINITY; jbe@0: box->lat_max = -INFINITY; jbe@0: box->lon_min = 0; jbe@0: box->lon_max = 0; jbe@0: } jbe@0: jbe@0: /* check if point is inside a box */ jbe@0: static bool pgl_point_in_box(pgl_point *point, pgl_box *box) { jbe@0: return ( jbe@0: point->lat >= box->lat_min && point->lat <= box->lat_max && ( jbe@0: (box->lon_min > box->lon_max) ? ( jbe@0: /* box crosses 180th meridian */ jbe@0: point->lon >= box->lon_min || point->lon <= box->lon_max jbe@0: ) : ( jbe@0: /* box does not cross the 180th meridian */ jbe@0: point->lon >= box->lon_min && point->lon <= box->lon_max jbe@0: ) jbe@0: ) jbe@0: ); jbe@0: } jbe@0: jbe@0: /* check if two boxes overlap */ jbe@0: static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) { jbe@0: return ( jbe@0: box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */ jbe@0: ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) && jbe@0: ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && ( jbe@0: ( jbe@0: /* check if one and only one box crosses the 180th meridian */ jbe@0: ((box1->lon_min > box1->lon_max) ? 1 : 0) ^ jbe@0: ((box2->lon_min > box2->lon_max) ? 1 : 0) jbe@0: ) ? ( jbe@0: /* exactly one box crosses the 180th meridian */ jbe@0: box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min || jbe@0: box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max jbe@0: ) : ( jbe@0: /* no box or both boxes cross the 180th meridian */ jbe@0: ( jbe@0: (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) && jbe@0: (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max) jbe@0: ) || jbe@0: /* handle W180 == E180 */ jbe@0: ( box1->lon_min == -180 && box2->lon_max == 180 ) || jbe@0: ( box2->lon_min == -180 && box1->lon_max == 180 ) jbe@0: ) jbe@0: ) jbe@0: ); jbe@0: } jbe@0: jbe@0: /* check unambiguousness of east/west orientation of cluster entries and set jbe@0: bounding circle of cluster */ jbe@0: static bool pgl_finalize_cluster(pgl_cluster *cluster) { jbe@0: int i, j; /* i: index of entry, j: index of point in entry */ jbe@0: int npoints; /* number of points in entry */ jbe@0: int total_npoints = 0; /* total number of points in cluster */ jbe@0: pgl_point *points; /* points in entry */ jbe@0: int lon_dir; /* first point of entry west (-1) or east (+1) */ jbe@0: double lon_break = 0; /* antipodal longitude of first point in entry */ jbe@0: double lon_min, lon_max; /* covered longitude range of entry */ jbe@0: double value; /* temporary variable */ jbe@0: /* reset bounding circle center to empty circle at 0/0 coordinates */ jbe@0: cluster->bounding.center.lat = 0; jbe@0: cluster->bounding.center.lon = 0; jbe@0: cluster->bounding.radius = -INFINITY; jbe@0: /* if cluster is not empty */ jbe@0: if (cluster->nentries != 0) { jbe@0: /* iterate over all cluster entries and ensure they each cover a longitude jbe@0: range less than 180 degrees */ jbe@0: for (i=0; inentries; i++) { jbe@0: /* get properties of entry */ jbe@0: npoints = cluster->entries[i].npoints; jbe@0: points = PGL_ENTRY_POINTS(cluster, i); jbe@0: /* get longitude of first point of entry */ jbe@0: value = points[0].lon; jbe@0: /* initialize lon_min and lon_max with longitude of first point */ jbe@0: lon_min = value; jbe@0: lon_max = value; jbe@0: /* determine east/west orientation of first point and calculate antipodal jbe@0: longitude (Note: rounding required here) */ jbe@0: if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); } jbe@0: else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); } jbe@0: else lon_dir = 0; jbe@0: /* iterate over all other points in entry */ jbe@0: for (j=1; jlon_break) value = pgl_round(value - 360); jbe@0: else if (lon_dir>0 && value lon_max) lon_max = value; jbe@0: /* return false if 180 degrees or more are covered */ jbe@0: if (lon_max - lon_min >= 180) return false; jbe@0: } jbe@0: } jbe@0: /* iterate over all points of all entries and calculate arbitrary center jbe@0: point for bounding circle (best if center point minimizes the radius, jbe@0: but some error is allowed here) */ jbe@0: for (i=0; inentries; i++) { jbe@0: /* get properties of entry */ jbe@0: npoints = cluster->entries[i].npoints; jbe@0: points = PGL_ENTRY_POINTS(cluster, i); jbe@0: /* check if first entry */ jbe@0: if (i==0) { jbe@0: /* get longitude of first point of first entry in whole cluster */ jbe@0: value = points[0].lon; jbe@0: /* initialize lon_min and lon_max with longitude of first point of jbe@0: first entry in whole cluster (used to determine if whole cluster jbe@0: covers a longitude range of 180 degrees or more) */ jbe@0: lon_min = value; jbe@0: lon_max = value; jbe@0: /* determine east/west orientation of first point and calculate jbe@0: antipodal longitude (Note: rounding not necessary here) */ jbe@0: if (value < 0) { lon_dir = -1; lon_break = value + 180; } jbe@0: else if (value > 0) { lon_dir = 1; lon_break = value - 180; } jbe@0: else lon_dir = 0; jbe@0: } jbe@0: /* iterate over all points in entry */ jbe@0: for (j=0; j lon_break) value -= 360; jbe@0: else if (lon_dir > 0 && value < lon_break) value += 360; jbe@0: if (value < lon_min) lon_min = value; jbe@0: else if (value > lon_max) lon_max = value; jbe@0: /* set bounding circle to cover whole earth if more than 180 degrees jbe@0: are covered */ jbe@0: if (lon_max - lon_min >= 180) { jbe@0: cluster->bounding.center.lat = 0; jbe@0: cluster->bounding.center.lon = 0; jbe@0: cluster->bounding.radius = INFINITY; jbe@0: return true; jbe@0: } jbe@0: /* add point to bounding circle center (for average calculation) */ jbe@0: cluster->bounding.center.lat += points[j].lat; jbe@0: cluster->bounding.center.lon += value; jbe@0: } jbe@0: /* count total number of points */ jbe@0: total_npoints += npoints; jbe@0: } jbe@0: /* determine average latitude and longitude of cluster */ jbe@0: cluster->bounding.center.lat /= total_npoints; jbe@0: cluster->bounding.center.lon /= total_npoints; jbe@0: /* normalize longitude of center of cluster bounding circle */ jbe@0: if (cluster->bounding.center.lon < -180) { jbe@0: cluster->bounding.center.lon += 360; jbe@0: } jbe@0: else if (cluster->bounding.center.lon > 180) { jbe@0: cluster->bounding.center.lon -= 360; jbe@0: } jbe@0: /* round bounding circle center (useful if it is used by other functions) */ jbe@0: cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat); jbe@0: cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon); jbe@0: /* calculate radius of bounding circle */ jbe@0: for (i=0; inentries; i++) { jbe@0: npoints = cluster->entries[i].npoints; jbe@0: points = PGL_ENTRY_POINTS(cluster, i); jbe@0: for (j=0; jbounding.center.lat, cluster->bounding.center.lon, jbe@0: points[j].lat, points[j].lon jbe@0: ); jbe@0: if (value > cluster->bounding.radius) cluster->bounding.radius = value; jbe@0: } jbe@0: } jbe@0: } jbe@0: /* return true (east/west orientation is unambiguous) */ jbe@0: return true; jbe@0: } jbe@0: jbe@0: /* check if point is inside cluster */ jbe@20: /* (if point is on perimeter, then true is returned if and only if jbe@20: strict == false) */ jbe@20: static bool pgl_point_in_cluster( jbe@20: pgl_point *point, jbe@20: pgl_cluster *cluster, jbe@20: bool strict jbe@20: ) { jbe@0: int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ jbe@0: int entrytype; /* type of entry */ jbe@0: int npoints; /* number of points in entry */ jbe@0: pgl_point *points; /* array of points in entry */ jbe@0: int lon_dir = 0; /* first vertex west (-1) or east (+1) */ jbe@0: double lon_break = 0; /* antipodal longitude of first vertex */ jbe@0: double lat0 = point->lat; /* latitude of point */ jbe@0: double lon0; /* (adjusted) longitude of point */ jbe@0: double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ jbe@0: double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ jbe@0: double lon; /* longitude of intersection */ jbe@0: int counter = 0; /* counter for intersections east of point */ jbe@0: /* iterate over all entries */ jbe@0: for (i=0; inentries; i++) { jbe@20: /* get type of entry */ jbe@0: entrytype = cluster->entries[i].entrytype; jbe@20: /* skip all entries but polygons if perimeters are excluded */ jbe@20: if (strict && entrytype != PGL_ENTRY_POLYGON) continue; jbe@20: /* get points of entry */ jbe@0: npoints = cluster->entries[i].npoints; jbe@0: points = PGL_ENTRY_POINTS(cluster, i); jbe@0: /* determine east/west orientation of first point of entry and calculate jbe@0: antipodal longitude */ jbe@0: lon_break = points[0].lon; jbe@0: if (lon_break < 0) { lon_dir = -1; lon_break += 180; } jbe@0: else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } jbe@0: else lon_dir = 0; jbe@0: /* get longitude of point */ jbe@0: lon0 = point->lon; jbe@0: /* consider longitude wrap-around for point */ jbe@0: if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360); jbe@0: else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360); jbe@0: /* iterate over all edges and vertices */ jbe@0: for (j=0; j lon_break) lon1 = pgl_round(lon1 - 360); jbe@16: else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360); jbe@16: } jbe@16: /* get latitude and longitude of next vertex */ jbe@0: lat2 = points[k].lat; jbe@0: lon2 = points[k].lon; jbe@16: /* consider longitude wrap-around for next vertex */ jbe@0: if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360); jbe@0: else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360); jbe@20: /* return if point is on horizontal (west to east) edge of polygon */ jbe@0: if ( jbe@0: lat0 == lat1 && lat0 == lat2 && jbe@0: ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) ) jbe@20: ) return !strict; jbe@0: /* check if edge crosses east/west line of point */ jbe@0: if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) { jbe@0: /* calculate longitude of intersection */ jbe@0: lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1); jbe@20: /* return if intersection goes (approximately) through point */ jbe@20: if (pgl_round(lon) == lon0) return !strict; jbe@0: /* count intersection if east of point and entry is polygon*/ jbe@0: if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++; jbe@0: } jbe@0: } jbe@0: } jbe@0: /* return true if number of intersections is odd */ jbe@0: return counter & 1; jbe@0: } jbe@0: jbe@20: /* check if all points of the second cluster are strictly inside the first jbe@20: cluster */ jbe@20: static inline bool pgl_all_cluster_points_strictly_in_cluster( jbe@16: pgl_cluster *outer, pgl_cluster *inner jbe@16: ) { jbe@16: int i, j; /* i: entry, j: point in entry */ jbe@16: int npoints; /* number of points in entry */ jbe@16: pgl_point *points; /* array of points in entry */ jbe@16: /* iterate over all entries of "inner" cluster */ jbe@16: for (i=0; inentries; i++) { jbe@16: /* get properties of entry */ jbe@16: npoints = inner->entries[i].npoints; jbe@16: points = PGL_ENTRY_POINTS(inner, i); jbe@16: /* iterate over all points in entry of "inner" cluster */ jbe@16: for (j=0; jnentries; i++) { jbe@16: /* get properties of entry */ jbe@16: npoints = inner->entries[i].npoints; jbe@16: points = PGL_ENTRY_POINTS(inner, i); jbe@16: /* iterate over all points in entry of "inner" cluster */ jbe@16: for (j=0; jnentries; i1++) { jbe@16: /* get properties of entry in cluster1 and skip points */ jbe@16: npoints1 = cluster1->entries[i1].npoints; jbe@16: if (npoints1 < 2) continue; jbe@16: entrytype1 = cluster1->entries[i1].entrytype; jbe@16: points1 = PGL_ENTRY_POINTS(cluster1, i1); jbe@16: /* determine east/west orientation of first point and calculate antipodal jbe@16: longitude */ jbe@16: lon_break1 = points1[0].lon; jbe@16: if (lon_break1 < 0) { jbe@16: lon_dir1 = -1; jbe@16: lon_break1 = pgl_round(lon_break1 + 180); jbe@16: } else if (lon_break1 > 0) { jbe@16: lon_dir1 = 1; jbe@16: lon_break1 = pgl_round(lon_break1 - 180); jbe@16: } else lon_dir1 = 0; jbe@16: /* iterate over all edges and vertices in cluster1 */ jbe@16: for (j1=0; j1lon_break1) lon11 = pgl_round(lon11-360); jbe@16: else if (lon_dir1>0 && lon11lon_break1) lon12 = pgl_round(lon12-360); jbe@16: else if (lon_dir1>0 && lon12nentries; i2++) { jbe@16: /* get points and number of points of entry in cluster2 */ jbe@16: npoints2 = cluster2->entries[i2].npoints; jbe@16: if (npoints2 < 2) continue; jbe@16: entrytype2 = cluster2->entries[i2].entrytype; jbe@16: points2 = PGL_ENTRY_POINTS(cluster2, i2); jbe@16: /* determine east/west orientation of first point and calculate antipodal jbe@16: longitude */ jbe@16: lon_break2 = points2[0].lon; jbe@16: if (lon_break2 < 0) { jbe@16: lon_dir2 = -1; jbe@16: lon_break2 = pgl_round(lon_break2 + 180); jbe@16: } else if (lon_break2 > 0) { jbe@16: lon_dir2 = 1; jbe@16: lon_break2 = pgl_round(lon_break2 - 180); jbe@16: } else lon_dir2 = 0; jbe@16: /* iterate over all edges and vertices in cluster2 */ jbe@16: for (j2=0; j2lon_break2) lon21 = pgl_round(lon21-360); jbe@16: else if (lon_dir2>0 && lon21lon_break2) lon22 = pgl_round(lon22-360); jbe@16: else if (lon_dir2>0 && lon22 360) { jbe@16: lon21 = pgl_round(lon21 - 360); jbe@16: lon22 = pgl_round(lon22 - 360); jbe@16: } else if (wrapvalue < -360) { jbe@16: lon21 = pgl_round(lon21 + 360); jbe@16: lon22 = pgl_round(lon22 + 360); jbe@16: } jbe@16: /* return true if segments overlap */ jbe@16: if ( jbe@16: pgl_lseg_crosses_line( jbe@16: lat11, lon11, lat12, lon12, jbe@20: lat21, lon21, lat22, lon22 jbe@16: ) && pgl_lseg_crosses_line( jbe@16: lat21, lon21, lat22, lon22, jbe@20: lat11, lon11, lat12, lon12 jbe@16: ) jbe@16: ) { jbe@16: return true; jbe@16: } jbe@16: } jbe@16: } jbe@16: } jbe@16: } jbe@16: /* otherwise return false */ jbe@16: return false; jbe@16: } jbe@16: jbe@16: /* check if second cluster is completely contained in first cluster */ jbe@16: static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) { jbe@20: if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false; jbe@20: if (pgl_any_cluster_points_in_cluster(inner, outer)) return false; jbe@20: if (pgl_outlines_overlap(outer, inner)) return false; jbe@16: return true; jbe@16: } jbe@16: jbe@16: /* check if two clusters overlap */ jbe@16: static bool pgl_clusters_overlap( jbe@16: pgl_cluster *cluster1, pgl_cluster *cluster2 jbe@16: ) { jbe@16: if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true; jbe@16: if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true; jbe@20: if (pgl_outlines_overlap(cluster1, cluster2)) return true; jbe@16: return false; jbe@16: } jbe@16: jbe@16: jbe@0: /* calculate (approximate) distance between point and cluster */ jbe@0: static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) { jbe@24: double comp; /* square of compression of meridians */ jbe@0: int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ jbe@0: int entrytype; /* type of entry */ jbe@0: int npoints; /* number of points in entry */ jbe@0: pgl_point *points; /* array of points in entry */ jbe@0: int lon_dir = 0; /* first vertex west (-1) or east (+1) */ jbe@0: double lon_break = 0; /* antipodal longitude of first vertex */ jbe@0: double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */ jbe@0: double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */ jbe@0: double lat0 = point->lat; /* latitude of point */ jbe@0: double lon0; /* (adjusted) longitude of point */ jbe@0: double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ jbe@0: double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ jbe@0: double s; /* scalar for vector calculations */ jbe@0: double dist; /* distance calculated in one step */ jbe@0: double min_dist = INFINITY; /* minimum distance */ jbe@0: /* distance is zero if point is contained in cluster */ jbe@20: if (pgl_point_in_cluster(point, cluster, false)) return 0; jbe@24: /* calculate (approximate) square compression of meridians */ jbe@24: /* TODO: use more exact formula based on WGS-84 */ jbe@24: comp = cos((lat0 / 180.0) * M_PI); jbe@24: comp *= comp; jbe@0: /* iterate over all entries */ jbe@0: for (i=0; inentries; i++) { jbe@0: /* get properties of entry */ jbe@0: entrytype = cluster->entries[i].entrytype; jbe@0: npoints = cluster->entries[i].npoints; jbe@0: points = PGL_ENTRY_POINTS(cluster, i); jbe@0: /* determine east/west orientation of first point of entry and calculate jbe@0: antipodal longitude */ jbe@0: lon_break = points[0].lon; jbe@0: if (lon_break < 0) { lon_dir = -1; lon_break += 180; } jbe@0: else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } jbe@0: else lon_dir = 0; jbe@0: /* determine covered longitude range */ jbe@0: for (j=0; j lon_break) lon1 -= 360; jbe@0: else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; jbe@0: /* update minimum and maximum longitude of polygon */ jbe@0: if (j == 0 || lon1 < lon_min) lon_min = lon1; jbe@0: if (j == 0 || lon1 > lon_max) lon_max = lon1; jbe@0: } jbe@0: /* adjust longitude wrap-around according to full longitude range */ jbe@0: lon_break = (lon_max + lon_min) / 2; jbe@0: if (lon_break < 0) { lon_dir = -1; lon_break += 180; } jbe@0: else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } jbe@0: /* get longitude of point */ jbe@0: lon0 = point->lon; jbe@0: /* consider longitude wrap-around for point */ jbe@0: if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360; jbe@0: else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360; jbe@0: /* iterate over all edges and vertices */ jbe@0: for (j=0; j lon_break) lon1 -= 360; jbe@16: else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; jbe@16: } jbe@0: /* calculate distance to vertex */ jbe@0: dist = pgl_distance(lat0, lon0, lat1, lon1); jbe@0: /* store calculated distance if smallest */ jbe@0: if (dist < min_dist) min_dist = dist; jbe@0: /* calculate index of next vertex */ jbe@0: k = (j+1) % npoints; jbe@0: /* skip last edge unless entry is (closed) outline or polygon */ jbe@0: if ( jbe@0: k == 0 && jbe@0: entrytype != PGL_ENTRY_OUTLINE && jbe@0: entrytype != PGL_ENTRY_POLYGON jbe@0: ) continue; jbe@16: /* get latitude and longitude of next vertex */ jbe@0: lat2 = points[k].lat; jbe@0: lon2 = points[k].lon; jbe@16: /* consider longitude wrap-around for next vertex */ jbe@0: if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360; jbe@0: else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360; jbe@0: /* go to next vertex and edge if edge is degenerated */ jbe@0: if (lat1 == lat2 && lon1 == lon2) continue; jbe@0: /* otherwise test if point can be projected onto edge of polygon */ jbe@0: s = ( jbe@24: ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) / jbe@24: ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1)) jbe@0: ); jbe@0: /* go to next vertex and edge if point cannot be projected */ jbe@0: if (!(s > 0 && s < 1)) continue; jbe@0: /* calculate distance from original point to projected point */ jbe@0: dist = pgl_distance( jbe@0: lat0, lon0, jbe@0: lat1 + s * (lat2-lat1), jbe@0: lon1 + s * (lon2-lon1) jbe@0: ); jbe@0: /* store calculated distance if smallest */ jbe@0: if (dist < min_dist) min_dist = dist; jbe@0: } jbe@0: } jbe@0: /* return minimum distance */ jbe@0: return min_dist; jbe@0: } jbe@0: jbe@16: /* calculate (approximate) distance between two clusters */ jbe@16: static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) { jbe@16: int i, j; /* i: entry, j: point in entry */ jbe@16: int npoints; /* number of points in entry */ jbe@16: pgl_point *points; /* array of points in entry */ jbe@16: double dist; /* distance calculated in one step */ jbe@16: double min_dist = INFINITY; /* minimum distance */ jbe@16: /* consider distance from each point in one cluster to the whole other */ jbe@16: for (i=0; inentries; i++) { jbe@16: npoints = cluster1->entries[i].npoints; jbe@16: points = PGL_ENTRY_POINTS(cluster1, i); jbe@16: for (j=0; jnentries; i++) { jbe@16: npoints = cluster2->entries[i].npoints; jbe@16: points = PGL_ENTRY_POINTS(cluster2, i); jbe@16: for (j=0; jlat_min > box->lat_max) return INFINITY; jbe@16: /* return zero if point is inside box */ jbe@0: if (pgl_point_in_box(point, box)) return 0; jbe@0: /* calculate delta longitude */ jbe@0: dlon = box->lon_max - box->lon_min; jbe@0: if (dlon < 0) dlon += 360; /* 180th meridian crossed */ jbe@16: /* if delta longitude is greater than 150 degrees, perform safe fall-back */ jbe@16: if (dlon > 150) return 0; jbe@16: /* calculate lower limit for distance (formula below requires dlon <= 150) */ jbe@16: /* TODO: provide better estimation function to improve performance */ jbe@16: distance = ( jbe@16: (1.0-1e-14) * pgl_distance( jbe@16: point->lat, jbe@16: point->lon, jbe@16: (box->lat_min + box->lat_max) / 2, jbe@16: box->lon_min + dlon/2 jbe@16: ) - pgl_distance( jbe@16: box->lat_min, box->lon_min, jbe@16: box->lat_max, box->lon_max jbe@16: ) jbe@16: ); jbe@16: /* truncate negative results to zero */ jbe@16: if (distance <= 0) distance = 0; jbe@16: /* return result */ jbe@16: return distance; jbe@0: } jbe@0: jbe@0: jbe@16: /*-------------------------------------------------* jbe@16: * geographic index based on space-filling curve * jbe@16: *-------------------------------------------------*/ jbe@0: jbe@0: /* number of bytes used for geographic (center) position in keys */ jbe@0: #define PGL_KEY_LATLON_BYTELEN 7 jbe@0: jbe@0: /* maximum reference value for logarithmic size of geographic objects */ jbe@0: #define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */ jbe@0: jbe@0: /* pointer to index key (either pgl_pointkey or pgl_areakey) */ jbe@0: typedef unsigned char *pgl_keyptr; jbe@0: jbe@0: /* index key for points (objects with zero area) on the spheroid */ jbe@0: /* bit 0..55: interspersed bits of latitude and longitude, jbe@0: bit 56..57: always zero, jbe@0: bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */ jbe@0: typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1]; jbe@0: jbe@0: /* index key for geographic objects on spheroid with area greater than zero */ jbe@0: /* bit 0..55: interspersed bits of latitude and longitude of center point, jbe@0: bit 56: always set to 1, jbe@0: bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.), jbe@0: bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to jbe@0: PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth jbe@0: = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL jbe@0: (with interspersed bits = 0 and node depth = 0) for keys which jbe@0: cover both empty and non-empty objects */ jbe@0: jbe@0: typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2]; jbe@0: jbe@0: /* helper macros for reading/writing index keys */ jbe@0: #define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN jbe@0: #define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1) jbe@0: #define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8) jbe@0: #define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1) jbe@0: #define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1) jbe@0: #define PGL_AREAKEY_TYPEMASK 0x80 jbe@0: #define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8))) jbe@0: #define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \ jbe@0: ( PGL_KEY_LATLONBIT(key1, n) ^ \ jbe@0: PGL_KEY_LATLONBIT(key2, n) ) jbe@0: #define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ jbe@0: PGL_AREAKEY_TYPEMASK) jbe@0: #define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ jbe@0: (PGL_AREAKEY_TYPEMASK-1)) jbe@0: #define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET]) jbe@0: #define PGL_KEY_OBJSIZE_EMPTY 126 jbe@0: #define PGL_KEY_OBJSIZE_UNIVERSAL 127 jbe@0: #define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \ jbe@0: (key)[PGL_KEY_OBJSIZE_OFFSET] == \ jbe@0: PGL_KEY_OBJSIZE_EMPTY ) jbe@0: #define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \ jbe@0: (key)[PGL_KEY_OBJSIZE_OFFSET] == \ jbe@0: PGL_KEY_OBJSIZE_UNIVERSAL ) jbe@0: jbe@0: /* set area key to match empty objects only */ jbe@0: static void pgl_key_set_empty(pgl_keyptr key) { jbe@0: memset(key, 0, sizeof(pgl_areakey)); jbe@0: /* Note: setting node depth to maximum is required for picksplit function */ jbe@0: key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; jbe@0: key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY; jbe@0: } jbe@0: jbe@0: /* set area key to match any object (including empty objects) */ jbe@0: static void pgl_key_set_universal(pgl_keyptr key) { jbe@0: memset(key, 0, sizeof(pgl_areakey)); jbe@0: key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK; jbe@0: key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL; jbe@0: } jbe@0: jbe@0: /* convert a point on earth into a max-depth key to be used in index */ jbe@0: static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) { jbe@0: double lat = point->lat; jbe@0: double lon = point->lon; jbe@0: int i; jbe@0: /* clear latitude and longitude bits */ jbe@0: memset(key, 0, PGL_KEY_LATLON_BYTELEN); jbe@0: /* set node depth to maximum and type bit to zero */ jbe@0: key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH; jbe@0: /* iterate over all latitude/longitude bit pairs */ jbe@0: for (i=0; i= 0) { jbe@0: key[i/4] |= 0x80 >> (2*(i%4)); jbe@0: lat *= 2; lat -= 90; jbe@0: } else { jbe@0: lat *= 2; lat += 90; jbe@0: } jbe@0: /* determine longitude bit */ jbe@0: if (lon >= 0) { jbe@0: key[i/4] |= 0x80 >> (2*(i%4)+1); jbe@0: lon *= 2; lon -= 180; jbe@0: } else { jbe@0: lon *= 2; lon += 180; jbe@0: } jbe@0: } jbe@0: } jbe@0: jbe@0: /* convert a circle on earth into a max-depth key to be used in an index */ jbe@0: static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) { jbe@0: /* handle special case of empty circle */ jbe@0: if (circle->radius < 0) { jbe@0: pgl_key_set_empty(key); jbe@0: return; jbe@0: } jbe@0: /* perform same action as for point keys */ jbe@0: pgl_point_to_key(&(circle->center), key); jbe@0: /* but overwrite type and node depth to fit area index key */ jbe@0: key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; jbe@0: /* check if radius is greater than (or equal to) reference size */ jbe@0: /* (treat equal values as greater values for numerical safety) */ jbe@0: if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) { jbe@0: /* if yes, set logarithmic size to zero */ jbe@0: key[PGL_KEY_OBJSIZE_OFFSET] = 0; jbe@0: } else { jbe@0: /* otherwise, determine logarithmic size iteratively */ jbe@0: /* (one step is equivalent to a factor of sqrt(2)) */ jbe@0: double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2; jbe@0: int objsize = 1; jbe@0: while (objsize < PGL_AREAKEY_MAXOBJSIZE) { jbe@0: /* stop when radius is greater than (or equal to) adjusted reference */ jbe@0: /* (treat equal values as greater values for numerical safety) */ jbe@0: if (circle->radius >= reference) break; jbe@0: reference /= M_SQRT2; jbe@0: objsize++; jbe@0: } jbe@0: /* set logarithmic size to determined value */ jbe@0: key[PGL_KEY_OBJSIZE_OFFSET] = objsize; jbe@0: } jbe@0: } jbe@0: jbe@0: /* check if one key is subkey of another key or vice versa */ jbe@0: static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) { jbe@0: int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ jbe@0: /* determine smallest depth */ jbe@0: int depth1 = PGL_KEY_NODEDEPTH(key1); jbe@0: int depth2 = PGL_KEY_NODEDEPTH(key2); jbe@0: int depth = (depth1 < depth2) ? depth1 : depth2; jbe@0: /* check if keys are area keys (assuming that both keys have same type) */ jbe@0: if (PGL_KEY_IS_AREAKEY(key1)) { jbe@0: int j = 0; /* bit offset for logarithmic object size bits */ jbe@0: int k = 0; /* bit offset for latitude and longitude */ jbe@0: /* fetch logarithmic object size information */ jbe@0: int objsize1 = PGL_KEY_OBJSIZE(key1); jbe@0: int objsize2 = PGL_KEY_OBJSIZE(key2); jbe@0: /* handle special cases for empty objects (universal and empty keys) */ jbe@0: if ( jbe@0: objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL || jbe@0: objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL jbe@0: ) return true; jbe@0: if ( jbe@0: objsize1 == PGL_KEY_OBJSIZE_EMPTY || jbe@0: objsize2 == PGL_KEY_OBJSIZE_EMPTY jbe@0: ) return objsize1 == objsize2; jbe@0: /* iterate through key bits */ jbe@0: for (i=0; i j) || jbe@0: (objsize2 <= j && objsize1 > j) jbe@0: ) { jbe@0: /* bit differs, therefore keys are in separate branches */ jbe@0: return false; jbe@0: } jbe@0: /* increase bit counter for object size bits */ jbe@0: j++; jbe@0: } jbe@0: /* all other bits describe latitude and longitude */ jbe@0: else { jbe@0: /* check if bit differs in both keys */ jbe@0: if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) { jbe@0: /* bit differs, therefore keys are in separate branches */ jbe@0: return false; jbe@0: } jbe@0: /* increase bit counter for latitude/longitude bits */ jbe@0: k++; jbe@0: } jbe@0: } jbe@0: } jbe@0: /* if not, keys are point keys */ jbe@0: else { jbe@0: /* iterate through key bits */ jbe@0: for (i=0; i PGL_AREAKEY_MAXOBJSIZE || jbe@0: objsize2 > PGL_AREAKEY_MAXOBJSIZE jbe@0: ) { jbe@0: if ( jbe@0: objsize1 == PGL_KEY_OBJSIZE_EMPTY && jbe@0: objsize2 == PGL_KEY_OBJSIZE_EMPTY jbe@0: ) pgl_key_set_empty(dst); jbe@0: else pgl_key_set_universal(dst); jbe@0: return; jbe@0: } jbe@0: /* iterate through key bits */ jbe@0: for (i=0; i= j && objsize2 >= j) { jbe@0: /* set objsize in destination buffer to indicate that size bit is jbe@0: unset in destination buffer at the current bit position */ jbe@0: dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j; jbe@0: } jbe@0: /* break if object size bit is set in one key only */ jbe@0: else if (objsize1 >= j || objsize2 >= j) break; jbe@0: } jbe@0: /* all other bits describe latitude and longitude */ jbe@0: else { jbe@0: /* break if bit differs in both keys */ jbe@0: if (PGL_KEY_LATLONBIT(dst, k)) { jbe@0: if (!PGL_KEY_LATLONBIT(src, k)) break; jbe@0: /* but set bit in destination buffer if bit is set in both keys */ jbe@0: dstbuf[k/8] |= 0x80 >> (k%8); jbe@0: } else if (PGL_KEY_LATLONBIT(src, k)) break; jbe@0: /* increase bit counter for latitude/longitude bits */ jbe@0: k++; jbe@0: } jbe@0: } jbe@0: /* set common node depth and type bit (type bit = 1) */ jbe@0: dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i; jbe@0: /* copy contents of destination buffer to first key */ jbe@0: memcpy(dst, dstbuf, sizeof(pgl_areakey)); jbe@0: } jbe@0: /* if not, keys are point keys */ jbe@0: else { jbe@0: pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */ jbe@0: /* iterate through key bits */ jbe@0: for (i=0; i> (i%8); jbe@0: } else if (PGL_KEY_LATLONBIT(src, i)) break; jbe@0: } jbe@0: /* set common node depth (type bit = 0) */ jbe@0: dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i; jbe@0: /* copy contents of destination buffer to first key */ jbe@0: memcpy(dst, dstbuf, sizeof(pgl_pointkey)); jbe@0: } jbe@0: } jbe@0: jbe@0: /* determine center(!) boundaries and radius estimation of index key */ jbe@0: static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) { jbe@0: int i; jbe@0: /* determine node depth */ jbe@0: int depth = PGL_KEY_NODEDEPTH(key); jbe@0: /* center point of possible result */ jbe@0: double lat = 0; jbe@0: double lon = 0; jbe@0: /* maximum distance of real center point from key center */ jbe@0: double dlat = 90; jbe@0: double dlon = 180; jbe@0: /* maximum radius of contained objects */ jbe@0: double radius = 0; /* always return zero for point index keys */ jbe@0: /* check if key is area key */ jbe@0: if (PGL_KEY_IS_AREAKEY(key)) { jbe@0: /* get logarithmic object size */ jbe@0: int objsize = PGL_KEY_OBJSIZE(key); jbe@0: /* handle special cases for empty objects (universal and empty keys) */ jbe@0: if (objsize == PGL_KEY_OBJSIZE_EMPTY) { jbe@0: pgl_box_set_empty(box); jbe@0: return 0; jbe@0: } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) { jbe@0: box->lat_min = -90; jbe@0: box->lat_max = 90; jbe@0: box->lon_min = -180; jbe@0: box->lon_max = 180; jbe@0: return 0; /* any value >= 0 would do */ jbe@0: } jbe@0: /* calculate maximum possible radius of objects covered by the given key */ jbe@0: if (objsize == 0) radius = INFINITY; jbe@0: else { jbe@0: radius = PGL_AREAKEY_REFOBJSIZE; jbe@0: while (--objsize) radius /= M_SQRT2; jbe@0: } jbe@0: /* iterate over latitude and longitude bits in key */ jbe@0: /* (every second bit is a latitude or longitude bit) */ jbe@0: for (i=0; ilat_min = lat - dlat; jbe@0: box->lat_max = lat + dlat; jbe@0: box->lon_min = lon - dlon; jbe@0: box->lon_max = lon + dlon; jbe@0: /* return radius (as a function return value) */ jbe@0: return radius; jbe@0: } jbe@0: jbe@0: /* estimator function for distance between point and index key */ jbe@16: /* always returns a smaller value than actually correct or zero */ jbe@0: static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) { jbe@0: pgl_box box; /* center(!) bounding box of area index key */ jbe@0: /* calculate center(!) bounding box and maximum radius of objects covered jbe@0: by area index key (radius is zero for point index keys) */ jbe@0: double distance = pgl_key_to_box(key, &box); jbe@0: /* calculate estimated distance between bounding box of center point of jbe@0: indexed object and point passed as second argument, then substract maximum jbe@0: radius of objects covered by index key */ jbe@16: distance = pgl_estimate_point_box_distance(point, &box) - distance; jbe@0: /* truncate negative results to zero */ jbe@0: if (distance <= 0) distance = 0; jbe@0: /* return result */ jbe@0: return distance; jbe@0: } jbe@0: jbe@0: jbe@0: /*---------------------------------* jbe@0: * helper functions for text I/O * jbe@0: *---------------------------------*/ jbe@0: jbe@0: #define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */ jbe@0: jbe@0: /* convert floating point number to string (round-trip safe) */ jbe@0: static void pgl_print_float(char *buf, double flt) { jbe@0: /* check if number is integral */ jbe@0: if (trunc(flt) == flt) { jbe@0: /* for integral floats use maximum precision */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); jbe@0: } else { jbe@0: /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt); jbe@0: if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt); jbe@0: if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); jbe@0: } jbe@0: } jbe@0: jbe@0: /* convert latitude floating point number (in degrees) to string */ jbe@0: static void pgl_print_lat(char *buf, double lat) { jbe@0: if (signbit(lat)) { jbe@0: /* treat negative latitudes (including -0) as south */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat); jbe@0: } else { jbe@0: /* treat positive latitudes (including +0) as north */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat); jbe@0: } jbe@0: } jbe@0: jbe@0: /* convert longitude floating point number (in degrees) to string */ jbe@0: static void pgl_print_lon(char *buf, double lon) { jbe@0: if (signbit(lon)) { jbe@0: /* treat negative longitudes (including -0) as west */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon); jbe@0: } else { jbe@0: /* treat positive longitudes (including +0) as east */ jbe@0: snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon); jbe@0: } jbe@0: } jbe@0: jbe@0: /* bit masks used as return value of pgl_scan() function */ jbe@0: #define PGL_SCAN_NONE 0 /* no value has been parsed */ jbe@0: #define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */ jbe@0: #define PGL_SCAN_LON (1<<1) /* longitude has been parsed */ jbe@0: #define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */ jbe@0: jbe@0: /* parse a coordinate (can be latitude or longitude) */ jbe@0: static int pgl_scan(char **str, double *lat, double *lon) { jbe@0: double val; jbe@0: int len; jbe@0: if ( jbe@0: sscanf(*str, " N %lf %n", &val, &len) || jbe@0: sscanf(*str, " n %lf %n", &val, &len) jbe@0: ) { jbe@0: *str += len; *lat = val; return PGL_SCAN_LAT; jbe@0: } jbe@0: if ( jbe@0: sscanf(*str, " S %lf %n", &val, &len) || jbe@0: sscanf(*str, " s %lf %n", &val, &len) jbe@0: ) { jbe@0: *str += len; *lat = -val; return PGL_SCAN_LAT; jbe@0: } jbe@0: if ( jbe@0: sscanf(*str, " E %lf %n", &val, &len) || jbe@0: sscanf(*str, " e %lf %n", &val, &len) jbe@0: ) { jbe@0: *str += len; *lon = val; return PGL_SCAN_LON; jbe@0: } jbe@0: if ( jbe@0: sscanf(*str, " W %lf %n", &val, &len) || jbe@0: sscanf(*str, " w %lf %n", &val, &len) jbe@0: ) { jbe@0: *str += len; *lon = -val; return PGL_SCAN_LON; jbe@0: } jbe@0: return PGL_SCAN_NONE; jbe@0: } jbe@0: jbe@0: jbe@0: /*-----------------* jbe@0: * SQL functions * jbe@0: *-----------------*/ jbe@0: jbe@0: /* Note: These function names use "epoint", "ebox", etc. notation here instead jbe@0: of "point", "box", etc. in order to distinguish them from any previously jbe@0: defined functions. */ jbe@0: jbe@0: /* function needed for dummy types and/or not implemented features */ jbe@0: PG_FUNCTION_INFO_V1(pgl_notimpl); jbe@0: Datum pgl_notimpl(PG_FUNCTION_ARGS) { jbe@0: ereport(ERROR, (errmsg("not implemented by pgLatLon"))); jbe@0: } jbe@0: jbe@0: /* set point to latitude and longitude (including checks) */ jbe@0: static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) { jbe@0: /* reject infinite or NaN values */ jbe@0: if (!isfinite(lat) || !isfinite(lon)) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_DATA_EXCEPTION), jbe@0: errmsg("epoint requires finite coordinates") jbe@0: )); jbe@0: } jbe@0: /* check latitude bounds */ jbe@0: if (lat < -90) { jbe@0: ereport(WARNING, (errmsg("latitude exceeds south pole"))); jbe@0: lat = -90; jbe@0: } else if (lat > 90) { jbe@0: ereport(WARNING, (errmsg("latitude exceeds north pole"))); jbe@0: lat = 90; jbe@0: } jbe@0: /* check longitude bounds */ jbe@0: if (lon < -180) { jbe@0: ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized"))); jbe@0: lon += 360 - trunc(lon / 360) * 360; jbe@0: } else if (lon > 180) { jbe@0: ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized"))); jbe@0: lon -= 360 + trunc(lon / 360) * 360; jbe@0: } jbe@0: /* store rounded latitude/longitude values for round-trip safety */ jbe@0: point->lat = pgl_round(lat); jbe@0: point->lon = pgl_round(lon); jbe@0: } jbe@0: jbe@0: /* create point ("epoint" in SQL) from latitude and longitude */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_epoint); jbe@0: Datum pgl_create_epoint(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); jbe@0: pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1)); jbe@0: PG_RETURN_POINTER(point); jbe@0: } jbe@0: jbe@0: /* parse point ("epoint" in SQL) */ jbe@0: /* format: '[NS] [EW]' */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_in); jbe@0: Datum pgl_epoint_in(PG_FUNCTION_ARGS) { jbe@0: char *str = PG_GETARG_CSTRING(0); /* input string */ jbe@0: char *strptr = str; /* current position within string */ jbe@0: int done = 0; /* bit mask storing if latitude or longitude was read */ jbe@0: double lat, lon; /* parsed values as double precision floats */ jbe@0: pgl_point *point; /* return value (to be palloc'ed) */ jbe@0: /* parse two floats (each latitude or longitude) separated by white-space */ jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: if (strptr != str && isspace(strptr[-1])) { jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: } jbe@0: /* require end of string, and latitude and longitude parsed successfully */ jbe@0: if (strptr[0] || done != PGL_SCAN_LATLON) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), jbe@0: errmsg("invalid input syntax for type epoint: \"%s\"", str) jbe@0: )); jbe@0: } jbe@0: /* allocate memory for result */ jbe@0: point = (pgl_point *)palloc(sizeof(pgl_point)); jbe@0: /* set latitude and longitude (and perform checks) */ jbe@0: pgl_epoint_set_latlon(point, lat, lon); jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(point); jbe@0: } jbe@0: jbe@0: /* create box ("ebox" in SQL) that is empty */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_empty_ebox); jbe@0: Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: pgl_box_set_empty(box); jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* set box to given boundaries (including checks) */ jbe@0: static void pgl_ebox_set_boundaries( jbe@0: pgl_box *box, jbe@0: double lat_min, double lat_max, double lon_min, double lon_max jbe@0: ) { jbe@0: /* if minimum latitude is greater than maximum latitude, return empty box */ jbe@0: if (lat_min > lat_max) { jbe@0: pgl_box_set_empty(box); jbe@0: return; jbe@0: } jbe@0: /* otherwise reject infinite or NaN values */ jbe@0: if ( jbe@0: !isfinite(lat_min) || !isfinite(lat_max) || jbe@0: !isfinite(lon_min) || !isfinite(lon_max) jbe@0: ) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_DATA_EXCEPTION), jbe@0: errmsg("ebox requires finite coordinates") jbe@0: )); jbe@0: } jbe@0: /* check latitude bounds */ jbe@0: if (lat_max < -90) { jbe@0: ereport(WARNING, (errmsg("northern latitude exceeds south pole"))); jbe@0: lat_max = -90; jbe@0: } else if (lat_max > 90) { jbe@0: ereport(WARNING, (errmsg("northern latitude exceeds north pole"))); jbe@0: lat_max = 90; jbe@0: } jbe@0: if (lat_min < -90) { jbe@0: ereport(WARNING, (errmsg("southern latitude exceeds south pole"))); jbe@0: lat_min = -90; jbe@0: } else if (lat_min > 90) { jbe@0: ereport(WARNING, (errmsg("southern latitude exceeds north pole"))); jbe@0: lat_min = 90; jbe@0: } jbe@0: /* check if all longitudes are included */ jbe@0: if (lon_max - lon_min >= 360) { jbe@0: if (lon_max - lon_min > 360) ereport(WARNING, ( jbe@0: errmsg("longitude coverage greater than 360 degrees") jbe@0: )); jbe@0: lon_min = -180; jbe@0: lon_max = 180; jbe@0: } else { jbe@0: /* normalize longitude bounds */ jbe@0: if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360; jbe@0: else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360; jbe@0: if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360; jbe@0: else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360; jbe@0: } jbe@0: /* store rounded latitude/longitude values for round-trip safety */ jbe@0: box->lat_min = pgl_round(lat_min); jbe@0: box->lat_max = pgl_round(lat_max); jbe@0: box->lon_min = pgl_round(lon_min); jbe@0: box->lon_max = pgl_round(lon_max); jbe@0: /* ensure that rounding does not change orientation */ jbe@0: if (lon_min > lon_max && box->lon_min == box->lon_max) { jbe@0: box->lon_min = -180; jbe@0: box->lon_max = 180; jbe@0: } jbe@0: } jbe@0: jbe@0: /* create box ("ebox" in SQL) from min/max latitude and min/max longitude */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_ebox); jbe@0: Datum pgl_create_ebox(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: pgl_ebox_set_boundaries( jbe@0: box, jbe@0: PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), jbe@0: PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3) jbe@0: ); jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* create box ("ebox" in SQL) from two points ("epoint"s) */ jbe@0: /* (can not be used to cover a longitude range of more than 120 degrees) */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints); jbe@0: Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); jbe@0: pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: double lat_min, lat_max, lon_min, lon_max; jbe@0: double dlon; /* longitude range (delta longitude) */ jbe@0: /* order latitude and longitude boundaries */ jbe@0: if (point2->lat < point1->lat) { jbe@0: lat_min = point2->lat; jbe@0: lat_max = point1->lat; jbe@0: } else { jbe@0: lat_min = point1->lat; jbe@0: lat_max = point2->lat; jbe@0: } jbe@0: if (point2->lon < point1->lon) { jbe@0: lon_min = point2->lon; jbe@0: lon_max = point1->lon; jbe@0: } else { jbe@0: lon_min = point1->lon; jbe@0: lon_max = point2->lon; jbe@0: } jbe@0: /* calculate longitude range (round to avoid floating point errors) */ jbe@0: dlon = pgl_round(lon_max - lon_min); jbe@0: /* determine east-west direction */ jbe@0: if (dlon >= 240) { jbe@0: /* assume that 180th meridian is crossed and swap min/max longitude */ jbe@0: double swap = lon_min; lon_min = lon_max; lon_max = swap; jbe@0: } else if (dlon > 120) { jbe@0: /* unclear orientation since delta longitude > 120 */ jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_DATA_EXCEPTION), jbe@0: errmsg("can not determine east/west orientation for ebox") jbe@0: )); jbe@0: } jbe@0: /* use boundaries to setup box (and perform checks) */ jbe@0: pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* parse box ("ebox" in SQL) */ jbe@0: /* format: '[NS] [EW] [NS] [EW]' jbe@0: or: '[NS] [NS] [EW] [EW]' */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_in); jbe@0: Datum pgl_ebox_in(PG_FUNCTION_ARGS) { jbe@0: char *str = PG_GETARG_CSTRING(0); /* input string */ jbe@0: char *str_lower; /* lower case version of input string */ jbe@0: char *strptr; /* current position within string */ jbe@0: int valid; /* number of valid chars */ jbe@0: int done; /* specifies if latitude or longitude was read */ jbe@0: double val; /* temporary variable */ jbe@0: int lat_count = 0; /* count of latitude values parsed */ jbe@0: int lon_count = 0; /* count of longitufde values parsed */ jbe@0: double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */ jbe@0: pgl_box *box; /* return value (to be palloc'ed) */ jbe@0: /* lowercase input */ jbe@0: str_lower = psprintf("%s", str); jbe@0: for (strptr=str_lower; *strptr; strptr++) { jbe@0: if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; jbe@0: } jbe@0: /* reset reading position to start of (lowercase) string */ jbe@0: strptr = str_lower; jbe@0: /* check if empty box */ jbe@0: valid = 0; jbe@0: sscanf(strptr, " empty %n", &valid); jbe@0: if (valid && strptr[valid] == 0) { jbe@0: /* allocate and return empty box */ jbe@0: box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: pgl_box_set_empty(box); jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: /* demand four blocks separated by whitespace */ jbe@0: valid = 0; jbe@0: sscanf(strptr, " %*s %*s %*s %*s %n", &valid); jbe@0: /* if four blocks separated by whitespace exist, parse those blocks */ jbe@0: if (strptr[valid] == 0) while (strptr[0]) { jbe@0: /* parse either latitude or longitude (whichever found in input string) */ jbe@0: done = pgl_scan(&strptr, &val, &val); jbe@0: /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */ jbe@0: if (done == PGL_SCAN_LAT) { jbe@0: if (!lat_count) lat_min = val; else lat_max = val; jbe@0: lat_count++; jbe@0: } else if (done == PGL_SCAN_LON) { jbe@0: if (!lon_count) lon_min = val; else lon_max = val; jbe@0: lon_count++; jbe@0: } else { jbe@0: break; jbe@0: } jbe@0: } jbe@0: /* require end of string, and two latitude and two longitude values */ jbe@0: if (strptr[0] || lat_count != 2 || lon_count != 2) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), jbe@0: errmsg("invalid input syntax for type ebox: \"%s\"", str) jbe@0: )); jbe@0: } jbe@0: /* free lower case string */ jbe@0: pfree(str_lower); jbe@0: /* order boundaries (maximum greater than minimum) */ jbe@0: if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; } jbe@0: if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; } jbe@0: /* allocate memory for result */ jbe@0: box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: /* set boundaries (and perform checks) */ jbe@0: pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* set circle to given latitude, longitude, and radius (including checks) */ jbe@0: static void pgl_ecircle_set_latlon_radius( jbe@0: pgl_circle *circle, double lat, double lon, double radius jbe@0: ) { jbe@0: /* set center point (including checks) */ jbe@0: pgl_epoint_set_latlon(&(circle->center), lat, lon); jbe@0: /* handle non-positive radius */ jbe@0: if (isnan(radius)) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_DATA_EXCEPTION), jbe@0: errmsg("invalid radius for ecircle") jbe@0: )); jbe@0: } jbe@0: if (radius == 0) radius = 0; /* avoids -0 */ jbe@0: else if (radius < 0) { jbe@0: if (isfinite(radius)) { jbe@0: ereport(NOTICE, (errmsg("negative radius converted to minus infinity"))); jbe@0: } jbe@0: radius = -INFINITY; jbe@0: } jbe@0: /* store radius (round-trip safety is ensured by pgl_print_float) */ jbe@0: circle->radius = radius; jbe@0: } jbe@0: jbe@0: /* create circle ("ecircle" in SQL) from latitude, longitude, and radius */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_ecircle); jbe@0: Datum pgl_create_ecircle(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); jbe@0: pgl_ecircle_set_latlon_radius( jbe@0: circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2) jbe@0: ); jbe@0: PG_RETURN_POINTER(circle); jbe@0: } jbe@0: jbe@0: /* create circle ("ecircle" in SQL) from point ("epoint"), and radius */ jbe@0: PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint); jbe@0: Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: double radius = PG_GETARG_FLOAT8(1); jbe@0: pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); jbe@0: /* set latitude, longitude, radius (and perform checks) */ jbe@0: pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius); jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(circle); jbe@0: } jbe@0: jbe@0: /* parse circle ("ecircle" in SQL) */ jbe@0: /* format: '[NS] [EW] ' */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_in); jbe@0: Datum pgl_ecircle_in(PG_FUNCTION_ARGS) { jbe@0: char *str = PG_GETARG_CSTRING(0); /* input string */ jbe@0: char *strptr = str; /* current position within string */ jbe@0: double lat, lon, radius; /* parsed values as double precision flaots */ jbe@0: int valid = 0; /* number of valid chars */ jbe@0: int done = 0; /* stores if latitude and/or longitude was read */ jbe@0: pgl_circle *circle; /* return value (to be palloc'ed) */ jbe@0: /* demand three blocks separated by whitespace */ jbe@0: sscanf(strptr, " %*s %*s %*s %n", &valid); jbe@0: /* if three blocks separated by whitespace exist, parse those blocks */ jbe@0: if (strptr[valid] == 0) { jbe@0: /* parse latitude and longitude */ jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: /* parse radius (while incrementing strptr by number of bytes parsed) */ jbe@0: valid = 0; jbe@0: if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid; jbe@0: } jbe@0: /* require end of string and both latitude and longitude being parsed */ jbe@0: if (strptr[0] || done != PGL_SCAN_LATLON) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), jbe@0: errmsg("invalid input syntax for type ecircle: \"%s\"", str) jbe@0: )); jbe@0: } jbe@0: /* allocate memory for result */ jbe@0: circle = (pgl_circle *)palloc(sizeof(pgl_circle)); jbe@0: /* set latitude, longitude, radius (and perform checks) */ jbe@0: pgl_ecircle_set_latlon_radius(circle, lat, lon, radius); jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(circle); jbe@0: } jbe@0: jbe@0: /* parse cluster ("ecluster" in SQL) */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecluster_in); jbe@0: Datum pgl_ecluster_in(PG_FUNCTION_ARGS) { jbe@0: int i; jbe@0: char *str = PG_GETARG_CSTRING(0); /* input string */ jbe@0: char *str_lower; /* lower case version of input string */ jbe@0: char *strptr; /* pointer to current reading position of input */ jbe@0: int npoints_total = 0; /* total number of points in cluster */ jbe@0: int nentries = 0; /* total number of entries */ jbe@0: pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */ jbe@0: int entries_buflen = 4; /* maximum number of elements in entries array */ jbe@0: int valid; /* number of valid chars processed */ jbe@0: double lat, lon; /* latitude and longitude of parsed point */ jbe@0: int entrytype; /* current entry type */ jbe@0: int npoints; /* number of points in current entry */ jbe@0: pgl_point *points; /* array of pgl_point for pgl_newentry */ jbe@0: int points_buflen; /* maximum number of elements in points array */ jbe@0: int done; /* return value of pgl_scan function */ jbe@0: pgl_cluster *cluster; /* created cluster */ jbe@0: /* lowercase input */ jbe@0: str_lower = psprintf("%s", str); jbe@0: for (strptr=str_lower; *strptr; strptr++) { jbe@0: if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; jbe@0: } jbe@0: /* reset reading position to start of (lowercase) string */ jbe@0: strptr = str_lower; jbe@0: /* allocate initial buffer for entries */ jbe@0: entries = palloc(entries_buflen * sizeof(pgl_newentry)); jbe@0: /* parse until end of string */ jbe@0: while (strptr[0]) { jbe@0: /* require previous white-space or closing parenthesis before next token */ jbe@0: if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') { jbe@0: goto pgl_ecluster_in_error; jbe@0: } jbe@0: /* ignore token "empty" */ jbe@0: valid = 0; sscanf(strptr, " empty %n", &valid); jbe@0: if (valid) { strptr += valid; continue; } jbe@0: /* test for "point" token */ jbe@0: valid = 0; sscanf(strptr, " point ( %n", &valid); jbe@0: if (valid) { jbe@0: strptr += valid; jbe@0: entrytype = PGL_ENTRY_POINT; jbe@0: goto pgl_ecluster_in_type_ok; jbe@0: } jbe@0: /* test for "path" token */ jbe@0: valid = 0; sscanf(strptr, " path ( %n", &valid); jbe@0: if (valid) { jbe@0: strptr += valid; jbe@0: entrytype = PGL_ENTRY_PATH; jbe@0: goto pgl_ecluster_in_type_ok; jbe@0: } jbe@0: /* test for "outline" token */ jbe@0: valid = 0; sscanf(strptr, " outline ( %n", &valid); jbe@0: if (valid) { jbe@0: strptr += valid; jbe@0: entrytype = PGL_ENTRY_OUTLINE; jbe@0: goto pgl_ecluster_in_type_ok; jbe@0: } jbe@0: /* test for "polygon" token */ jbe@0: valid = 0; sscanf(strptr, " polygon ( %n", &valid); jbe@0: if (valid) { jbe@0: strptr += valid; jbe@0: entrytype = PGL_ENTRY_POLYGON; jbe@0: goto pgl_ecluster_in_type_ok; jbe@0: } jbe@0: /* error if no valid token found */ jbe@0: goto pgl_ecluster_in_error; jbe@0: pgl_ecluster_in_type_ok: jbe@0: /* check if pgl_newentry array needs to grow */ jbe@0: if (nentries == entries_buflen) { jbe@0: pgl_newentry *newbuf; jbe@0: entries_buflen *= 2; jbe@0: newbuf = palloc(entries_buflen * sizeof(pgl_newentry)); jbe@0: memcpy(newbuf, entries, nentries * sizeof(pgl_newentry)); jbe@0: pfree(entries); jbe@0: entries = newbuf; jbe@0: } jbe@0: /* reset number of points for current entry */ jbe@0: npoints = 0; jbe@0: /* allocate array for points */ jbe@0: points_buflen = 4; jbe@0: points = palloc(points_buflen * sizeof(pgl_point)); jbe@0: /* parse until closing parenthesis */ jbe@0: while (strptr[0] != ')') { jbe@0: /* error on unexpected end of string */ jbe@0: if (strptr[0] == 0) goto pgl_ecluster_in_error; jbe@0: /* mark neither latitude nor longitude as read */ jbe@0: done = PGL_SCAN_NONE; jbe@0: /* require white-space before second, third, etc. point */ jbe@0: if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error; jbe@0: /* scan latitude (or longitude) */ jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: /* require white-space before second coordinate */ jbe@0: if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error; jbe@0: /* scan longitude (or latitude) */ jbe@0: done |= pgl_scan(&strptr, &lat, &lon); jbe@0: /* error unless both latitude and longitude were parsed */ jbe@0: if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error; jbe@0: /* throw error if number of points is too high */ jbe@0: if (npoints_total == PGL_CLUSTER_MAXPOINTS) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), jbe@0: errmsg( jbe@0: "too many points for ecluster entry (maximum %i)", jbe@0: PGL_CLUSTER_MAXPOINTS jbe@0: ) jbe@0: )); jbe@0: } jbe@0: /* check if pgl_point array needs to grow */ jbe@0: if (npoints == points_buflen) { jbe@0: pgl_point *newbuf; jbe@0: points_buflen *= 2; jbe@0: newbuf = palloc(points_buflen * sizeof(pgl_point)); jbe@0: memcpy(newbuf, points, npoints * sizeof(pgl_point)); jbe@0: pfree(points); jbe@0: points = newbuf; jbe@0: } jbe@0: /* append point to pgl_point array (includes checks) */ jbe@0: pgl_epoint_set_latlon(&(points[npoints++]), lat, lon); jbe@0: /* increase total number of points */ jbe@0: npoints_total++; jbe@0: } jbe@0: /* error if entry has no points */ jbe@0: if (!npoints) goto pgl_ecluster_in_error; jbe@0: /* entries with one point are automatically of type "point" */ jbe@0: if (npoints == 1) entrytype = PGL_ENTRY_POINT; jbe@0: /* if entries have more than one point */ jbe@0: else { jbe@0: /* throw error if entry type is "point" */ jbe@0: if (entrytype == PGL_ENTRY_POINT) { jbe@0: ereport(ERROR, ( jbe@0: errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), jbe@0: errmsg("invalid input syntax for type ecluster (point entry with more than one point)") jbe@0: )); jbe@0: } jbe@0: /* coerce outlines and polygons with more than 2 points to be a path */ jbe@0: if (npoints == 2) entrytype = PGL_ENTRY_PATH; jbe@0: } jbe@0: /* append entry to pgl_newentry array */ jbe@0: entries[nentries].entrytype = entrytype; jbe@0: entries[nentries].npoints = npoints; jbe@0: entries[nentries].points = points; jbe@0: nentries++; jbe@0: /* consume closing parenthesis */ jbe@0: strptr++; jbe@0: /* consume white-space */ jbe@0: while (isspace(strptr[0])) strptr++; jbe@0: } jbe@0: /* free lower case string */ jbe@0: pfree(str_lower); jbe@0: /* create cluster from pgl_newentry array */ jbe@0: cluster = pgl_new_cluster(nentries, entries); jbe@0: /* free pgl_newentry array */ jbe@0: for (i=0; ilat); jbe@0: pgl_print_lon(lonstr, point->lon); jbe@0: PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr)); jbe@0: } jbe@0: jbe@0: /* convert box ("ebox") to string representation */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_out); jbe@0: Datum pgl_ebox_out(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); jbe@0: double lon_min = box->lon_min; jbe@0: double lon_max = box->lon_max; jbe@0: char lat_min_str[PGL_NUMBUFLEN]; jbe@0: char lat_max_str[PGL_NUMBUFLEN]; jbe@0: char lon_min_str[PGL_NUMBUFLEN]; jbe@0: char lon_max_str[PGL_NUMBUFLEN]; jbe@0: /* return string "empty" if box is set to be empty */ jbe@0: if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty"); jbe@0: /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */ jbe@0: /* (required since pgl_box_in orders the longitude boundaries) */ jbe@0: if (lon_min > lon_max) { jbe@0: if (lon_min + lon_max >= 0) lon_min -= 360; jbe@0: else lon_max += 360; jbe@0: } jbe@0: /* format and return result */ jbe@0: pgl_print_lat(lat_min_str, box->lat_min); jbe@0: pgl_print_lat(lat_max_str, box->lat_max); jbe@0: pgl_print_lon(lon_min_str, lon_min); jbe@0: pgl_print_lon(lon_max_str, lon_max); jbe@0: PG_RETURN_CSTRING(psprintf( jbe@0: "%s %s %s %s", jbe@0: lat_min_str, lon_min_str, lat_max_str, lon_max_str jbe@0: )); jbe@0: } jbe@0: jbe@0: /* convert circle ("ecircle") to string representation */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_out); jbe@0: Datum pgl_ecircle_out(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: char latstr[PGL_NUMBUFLEN]; jbe@0: char lonstr[PGL_NUMBUFLEN]; jbe@0: char radstr[PGL_NUMBUFLEN]; jbe@0: pgl_print_lat(latstr, circle->center.lat); jbe@0: pgl_print_lon(lonstr, circle->center.lon); jbe@0: pgl_print_float(radstr, circle->radius); jbe@0: PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr)); jbe@0: } jbe@0: jbe@0: /* convert cluster ("ecluster") to string representation */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecluster_out); jbe@0: Datum pgl_ecluster_out(PG_FUNCTION_ARGS) { jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); jbe@0: char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */ jbe@0: char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */ jbe@0: char ***strings; /* array of array of strings */ jbe@0: char *string; /* string of current token */ jbe@0: char *res, *resptr; /* result and pointer to current write position */ jbe@0: size_t reslen = 1; /* length of result (init with 1 for terminator) */ jbe@0: int npoints; /* number of points of current entry */ jbe@0: int i, j; /* i: entry, j: point in entry */ jbe@0: /* handle empty clusters */ jbe@0: if (cluster->nentries == 0) { jbe@0: /* free detoasted cluster (if copy) */ jbe@0: PG_FREE_IF_COPY(cluster, 0); jbe@0: /* return static result */ jbe@0: PG_RETURN_CSTRING("empty"); jbe@0: } jbe@0: /* allocate array of array of strings */ jbe@0: strings = palloc(cluster->nentries * sizeof(char **)); jbe@0: /* iterate over all entries in cluster */ jbe@0: for (i=0; inentries; i++) { jbe@0: /* get number of points in entry */ jbe@0: npoints = cluster->entries[i].npoints; jbe@0: /* allocate array of strings (one string for each point plus two extra) */ jbe@0: strings[i] = palloc((2 + npoints) * sizeof(char *)); jbe@0: /* determine opening string */ jbe@0: switch (cluster->entries[i].entrytype) { jbe@0: case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break; jbe@0: case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break; jbe@0: case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break; jbe@0: case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break; jbe@0: default: string = (i==0)?"unknown" :" unknown"; jbe@0: } jbe@0: /* use opening string as first string in array */ jbe@0: strings[i][0] = string; jbe@0: /* update result length (for allocating result string later) */ jbe@0: reslen += strlen(string); jbe@0: /* iterate over all points */ jbe@0: for (j=0; jnentries; i++) { jbe@0: npoints = cluster->entries[i].npoints; jbe@0: for (j=0; jlat = pq_getmsgfloat8(buf); jbe@0: point->lon = pq_getmsgfloat8(buf); jbe@0: PG_RETURN_POINTER(point); jbe@0: } jbe@0: jbe@0: /* binary input function for box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_recv); jbe@0: Datum pgl_ebox_recv(PG_FUNCTION_ARGS) { jbe@0: StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); jbe@0: pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); jbe@0: box->lat_min = pq_getmsgfloat8(buf); jbe@0: box->lat_max = pq_getmsgfloat8(buf); jbe@0: box->lon_min = pq_getmsgfloat8(buf); jbe@0: box->lon_max = pq_getmsgfloat8(buf); jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* binary input function for circle ("ecircle") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_recv); jbe@0: Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) { jbe@0: StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); jbe@0: circle->center.lat = pq_getmsgfloat8(buf); jbe@0: circle->center.lon = pq_getmsgfloat8(buf); jbe@0: circle->radius = pq_getmsgfloat8(buf); jbe@0: PG_RETURN_POINTER(circle); jbe@0: } jbe@0: jbe@0: /* TODO: binary receive function for cluster */ jbe@0: jbe@0: /* binary output function for point ("epoint") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_send); jbe@0: Datum pgl_epoint_send(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: StringInfoData buf; jbe@0: pq_begintypsend(&buf); jbe@0: pq_sendfloat8(&buf, point->lat); jbe@0: pq_sendfloat8(&buf, point->lon); jbe@0: PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); jbe@0: } jbe@0: jbe@0: /* binary output function for box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_send); jbe@0: Datum pgl_ebox_send(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); jbe@0: StringInfoData buf; jbe@0: pq_begintypsend(&buf); jbe@0: pq_sendfloat8(&buf, box->lat_min); jbe@0: pq_sendfloat8(&buf, box->lat_max); jbe@0: pq_sendfloat8(&buf, box->lon_min); jbe@0: pq_sendfloat8(&buf, box->lon_max); jbe@0: PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); jbe@0: } jbe@0: jbe@0: /* binary output function for circle ("ecircle") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_send); jbe@0: Datum pgl_ecircle_send(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: StringInfoData buf; jbe@0: pq_begintypsend(&buf); jbe@0: pq_sendfloat8(&buf, circle->center.lat); jbe@0: pq_sendfloat8(&buf, circle->center.lon); jbe@0: pq_sendfloat8(&buf, circle->radius); jbe@0: PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); jbe@0: } jbe@0: jbe@0: /* TODO: binary send functions for cluster */ jbe@0: jbe@0: /* cast point ("epoint") to box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox); jbe@0: Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_box *box = palloc(sizeof(pgl_box)); jbe@0: box->lat_min = point->lat; jbe@0: box->lat_max = point->lat; jbe@0: box->lon_min = point->lon; jbe@0: box->lon_max = point->lon; jbe@0: PG_RETURN_POINTER(box); jbe@0: } jbe@0: jbe@0: /* cast point ("epoint") to circle ("ecircle") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle); jbe@0: Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle = palloc(sizeof(pgl_box)); jbe@0: circle->center = *point; jbe@0: circle->radius = 0; jbe@0: PG_RETURN_POINTER(circle); jbe@0: } jbe@0: jbe@0: /* cast point ("epoint") to cluster ("ecluster") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster); jbe@0: Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_newentry entry; jbe@0: entry.entrytype = PGL_ENTRY_POINT; jbe@0: entry.npoints = 1; jbe@0: entry.points = point; jbe@0: PG_RETURN_POINTER(pgl_new_cluster(1, &entry)); jbe@0: } jbe@0: jbe@0: /* cast box ("ebox") to cluster ("ecluster") */ jbe@0: #define pgl_ebox_to_ecluster_macro(i, a, b) \ jbe@0: entries[i].entrytype = PGL_ENTRY_POLYGON; \ jbe@0: entries[i].npoints = 4; \ jbe@0: entries[i].points = points[i]; \ jbe@0: points[i][0].lat = box->lat_min; \ jbe@0: points[i][0].lon = (a); \ jbe@0: points[i][1].lat = box->lat_min; \ jbe@0: points[i][1].lon = (b); \ jbe@0: points[i][2].lat = box->lat_max; \ jbe@0: points[i][2].lon = (b); \ jbe@0: points[i][3].lat = box->lat_max; \ jbe@0: points[i][3].lon = (a); jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster); jbe@0: Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); jbe@0: double lon, dlon; jbe@0: int nentries; jbe@0: pgl_newentry entries[3]; jbe@0: pgl_point points[3][4]; jbe@0: if (box->lat_min > box->lat_max) { jbe@0: nentries = 0; jbe@0: } else if (box->lon_min > box->lon_max) { jbe@0: if (box->lon_min < 0) { jbe@0: lon = pgl_round((box->lon_min + 180) / 2.0); jbe@0: nentries = 3; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); jbe@0: pgl_ebox_to_ecluster_macro(1, lon, 180); jbe@0: pgl_ebox_to_ecluster_macro(2, -180, box->lon_max); jbe@0: } else if (box->lon_max > 0) { jbe@0: lon = pgl_round((box->lon_max - 180) / 2.0); jbe@0: nentries = 3; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); jbe@0: pgl_ebox_to_ecluster_macro(1, -180, lon); jbe@0: pgl_ebox_to_ecluster_macro(2, lon, box->lon_max); jbe@0: } else { jbe@0: nentries = 2; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); jbe@0: pgl_ebox_to_ecluster_macro(1, -180, box->lon_max); jbe@0: } jbe@0: } else { jbe@0: dlon = pgl_round(box->lon_max - box->lon_min); jbe@0: if (dlon < 180) { jbe@0: nentries = 1; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max); jbe@0: } else { jbe@0: lon = pgl_round((box->lon_min + box->lon_max) / 2.0); jbe@0: if ( jbe@0: pgl_round(lon - box->lon_min) < 180 && jbe@0: pgl_round(box->lon_max - lon) < 180 jbe@0: ) { jbe@0: nentries = 2; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); jbe@0: pgl_ebox_to_ecluster_macro(1, lon, box->lon_max); jbe@0: } else { jbe@0: nentries = 3; jbe@0: pgl_ebox_to_ecluster_macro(0, box->lon_min, -60); jbe@0: pgl_ebox_to_ecluster_macro(1, -60, 60); jbe@0: pgl_ebox_to_ecluster_macro(2, 60, box->lon_max); jbe@0: } jbe@0: } jbe@0: } jbe@0: PG_RETURN_POINTER(pgl_new_cluster(nentries, entries)); jbe@0: } jbe@0: jbe@0: /* extract latitude from point ("epoint") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_lat); jbe@0: Datum pgl_epoint_lat(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat); jbe@0: } jbe@0: jbe@0: /* extract longitude from point ("epoint") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_lon); jbe@0: Datum pgl_epoint_lon(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon); jbe@0: } jbe@0: jbe@0: /* extract minimum latitude from box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_lat_min); jbe@0: Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min); jbe@0: } jbe@0: jbe@0: /* extract maximum latitude from box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_lat_max); jbe@0: Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max); jbe@0: } jbe@0: jbe@0: /* extract minimum longitude from box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_lon_min); jbe@0: Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min); jbe@0: } jbe@0: jbe@0: /* extract maximum longitude from box ("ebox") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_lon_max); jbe@0: Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max); jbe@0: } jbe@0: jbe@0: /* extract center point from circle ("ecircle") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_center); jbe@0: Datum pgl_ecircle_center(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center)); jbe@0: } jbe@0: jbe@0: /* extract radius from circle ("ecircle") */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_radius); jbe@0: Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) { jbe@0: PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius); jbe@0: } jbe@0: jbe@0: /* check if point is inside box (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap); jbe@0: Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1); jbe@0: PG_RETURN_BOOL(pgl_point_in_box(point, box)); jbe@0: } jbe@0: jbe@0: /* check if point is inside circle (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap); jbe@0: Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: PG_RETURN_BOOL( jbe@0: pgl_distance( jbe@0: point->lat, point->lon, jbe@0: circle->center.lat, circle->center.lon jbe@0: ) <= circle->radius jbe@0: ); jbe@0: } jbe@0: jbe@0: /* check if point is inside cluster (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap); jbe@0: Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@16: bool retval; jbe@16: /* points outside bounding circle are always assumed to be non-overlapping jbe@16: (necessary for consistent table and index scans) */ jbe@16: if ( jbe@16: pgl_distance( jbe@16: point->lat, point->lon, jbe@16: cluster->bounding.center.lat, cluster->bounding.center.lon jbe@16: ) > cluster->bounding.radius jbe@16: ) retval = false; jbe@20: else retval = pgl_point_in_cluster(point, cluster, false); jbe@0: PG_FREE_IF_COPY(cluster, 1); jbe@0: PG_RETURN_BOOL(retval); jbe@0: } jbe@0: jbe@10: /* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */ jbe@10: PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap); jbe@10: Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) { jbe@10: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@10: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@10: bool retval = pgl_distance( jbe@10: point->lat, point->lon, jbe@10: cluster->bounding.center.lat, cluster->bounding.center.lon jbe@10: ) <= cluster->bounding.radius; jbe@10: PG_FREE_IF_COPY(cluster, 1); jbe@10: PG_RETURN_BOOL(retval); jbe@10: } jbe@10: jbe@0: /* check if two boxes overlap (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ebox_overlap); jbe@0: Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0); jbe@0: pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1); jbe@0: PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2)); jbe@0: } jbe@0: jbe@10: /* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */ jbe@10: PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap); jbe@10: Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) { jbe@10: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); jbe@10: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); jbe@10: PG_RETURN_BOOL( jbe@10: pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius jbe@10: ); jbe@10: } jbe@10: jbe@10: /* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */ jbe@10: PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap); jbe@10: Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) { jbe@10: pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); jbe@10: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@10: bool retval = pgl_estimate_point_box_distance( jbe@10: &cluster->bounding.center, jbe@10: box jbe@10: ) <= cluster->bounding.radius; jbe@10: PG_FREE_IF_COPY(cluster, 1); jbe@10: PG_RETURN_BOOL(retval); jbe@10: } jbe@10: jbe@0: /* check if two circles overlap (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_overlap); jbe@0: Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: PG_RETURN_BOOL( jbe@0: pgl_distance( jbe@0: circle1->center.lat, circle1->center.lon, jbe@0: circle2->center.lat, circle2->center.lon jbe@0: ) <= circle1->radius + circle2->radius jbe@0: ); jbe@0: } jbe@0: jbe@0: /* check if circle and cluster overlap (overlap operator "&&") in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap); jbe@0: Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@0: bool retval = ( jbe@0: pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius jbe@0: ); jbe@0: PG_FREE_IF_COPY(cluster, 1); jbe@0: PG_RETURN_BOOL(retval); jbe@0: } jbe@0: jbe@17: /* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */ jbe@10: PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap); jbe@10: Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) { jbe@10: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); jbe@10: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@10: bool retval = pgl_distance( jbe@10: circle->center.lat, circle->center.lon, jbe@10: cluster->bounding.center.lat, cluster->bounding.center.lon jbe@10: ) <= circle->radius + cluster->bounding.radius; jbe@10: PG_FREE_IF_COPY(cluster, 1); jbe@10: PG_RETURN_BOOL(retval); jbe@10: } jbe@10: jbe@16: /* check if two clusters overlap (overlap operator "&&") in SQL */ jbe@16: PG_FUNCTION_INFO_V1(pgl_ecluster_overlap); jbe@16: Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) { jbe@16: pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); jbe@16: pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@16: bool retval; jbe@16: /* clusters with non-touching bounding circles are always assumed to be jbe@16: non-overlapping (improves performance and is necessary for consistent jbe@16: table and index scans) */ jbe@16: if ( jbe@16: pgl_distance( jbe@16: cluster1->bounding.center.lat, cluster1->bounding.center.lon, jbe@16: cluster2->bounding.center.lat, cluster2->bounding.center.lon jbe@16: ) > cluster1->bounding.radius + cluster2->bounding.radius jbe@16: ) retval = false; jbe@16: else retval = pgl_clusters_overlap(cluster1, cluster2); jbe@16: PG_FREE_IF_COPY(cluster1, 0); jbe@16: PG_FREE_IF_COPY(cluster2, 1); jbe@16: PG_RETURN_BOOL(retval); jbe@16: } jbe@16: jbe@10: /* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */ jbe@10: PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap); jbe@10: Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) { jbe@10: pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); jbe@10: pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@10: bool retval = pgl_distance( jbe@10: cluster1->bounding.center.lat, cluster1->bounding.center.lon, jbe@10: cluster2->bounding.center.lat, cluster2->bounding.center.lon jbe@10: ) <= cluster1->bounding.radius + cluster2->bounding.radius; jbe@10: PG_FREE_IF_COPY(cluster1, 0); jbe@10: PG_FREE_IF_COPY(cluster2, 1); jbe@10: PG_RETURN_BOOL(retval); jbe@10: } jbe@10: jbe@16: /* check if second cluster is in first cluster (cont. operator "@>) in SQL */ jbe@16: PG_FUNCTION_INFO_V1(pgl_ecluster_contains); jbe@16: Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) { jbe@16: pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); jbe@16: pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@16: bool retval; jbe@16: /* clusters with non-touching bounding circles are always assumed to be jbe@16: non-overlapping (improves performance and is necessary for consistent jbe@16: table and index scans) */ jbe@16: if ( jbe@16: pgl_distance( jbe@16: outer->bounding.center.lat, outer->bounding.center.lon, jbe@16: inner->bounding.center.lat, inner->bounding.center.lon jbe@16: ) > outer->bounding.radius + inner->bounding.radius jbe@16: ) retval = false; jbe@16: else retval = pgl_cluster_in_cluster(outer, inner); jbe@16: PG_FREE_IF_COPY(outer, 0); jbe@16: PG_FREE_IF_COPY(inner, 1); jbe@16: PG_RETURN_BOOL(retval); jbe@16: } jbe@16: jbe@0: /* calculate distance between two points ("<->" operator) in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_distance); jbe@0: Datum pgl_epoint_distance(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); jbe@0: PG_RETURN_FLOAT8(pgl_distance( jbe@0: point1->lat, point1->lon, point2->lat, point2->lon jbe@0: )); jbe@0: } jbe@0: jbe@0: /* calculate point to circle distance ("<->" operator) in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance); jbe@0: Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: double distance = pgl_distance( jbe@0: point->lat, point->lon, circle->center.lat, circle->center.lon jbe@0: ) - circle->radius; jbe@0: PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); jbe@0: } jbe@0: jbe@0: /* calculate point to cluster distance ("<->" operator) in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance); jbe@0: Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) { jbe@0: pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@0: double distance = pgl_point_cluster_distance(point, cluster); jbe@0: PG_FREE_IF_COPY(cluster, 1); jbe@0: PG_RETURN_FLOAT8(distance); jbe@0: } jbe@0: jbe@0: /* calculate distance between two circles ("<->" operator) in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_distance); jbe@0: Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: double distance = pgl_distance( jbe@0: circle1->center.lat, circle1->center.lon, jbe@0: circle2->center.lat, circle2->center.lon jbe@0: ) - (circle1->radius + circle2->radius); jbe@0: PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); jbe@0: } jbe@0: jbe@0: /* calculate circle to cluster distance ("<->" operator) in SQL */ jbe@0: PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance); jbe@0: Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) { jbe@0: pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@0: double distance = ( jbe@0: pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius jbe@0: ); jbe@0: PG_FREE_IF_COPY(cluster, 1); jbe@0: PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); jbe@0: } jbe@0: jbe@16: /* calculate distance between two clusters ("<->" operator) in SQL */ jbe@16: PG_FUNCTION_INFO_V1(pgl_ecluster_distance); jbe@16: Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) { jbe@16: pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); jbe@16: pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@16: double retval = pgl_cluster_distance(cluster1, cluster2); jbe@16: PG_FREE_IF_COPY(cluster1, 0); jbe@16: PG_FREE_IF_COPY(cluster2, 1); jbe@16: PG_RETURN_FLOAT8(retval); jbe@16: } jbe@16: jbe@0: jbe@0: /*-----------------------------------------------------------* jbe@0: * B-tree comparison operators and index support functions * jbe@0: *-----------------------------------------------------------*/ jbe@0: jbe@0: /* macro for a B-tree operator (without detoasting) */ jbe@0: #define PGL_BTREE_OPER(func, type, cmpfunc, oper) \ jbe@0: PG_FUNCTION_INFO_V1(func); \ jbe@0: Datum func(PG_FUNCTION_ARGS) { \ jbe@0: type *a = (type *)PG_GETARG_POINTER(0); \ jbe@0: type *b = (type *)PG_GETARG_POINTER(1); \ jbe@0: PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \ jbe@0: } jbe@0: jbe@0: /* macro for a B-tree comparison function (without detoasting) */ jbe@0: #define PGL_BTREE_CMP(func, type, cmpfunc) \ jbe@0: PG_FUNCTION_INFO_V1(func); \ jbe@0: Datum func(PG_FUNCTION_ARGS) { \ jbe@0: type *a = (type *)PG_GETARG_POINTER(0); \ jbe@0: type *b = (type *)PG_GETARG_POINTER(1); \ jbe@0: PG_RETURN_INT32(cmpfunc(a, b)); \ jbe@0: } jbe@0: jbe@0: /* macro for a B-tree operator (with detoasting) */ jbe@0: #define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \ jbe@0: PG_FUNCTION_INFO_V1(func); \ jbe@0: Datum func(PG_FUNCTION_ARGS) { \ jbe@0: bool res; \ jbe@0: type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ jbe@0: type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ jbe@0: res = cmpfunc(a, b) oper 0; \ jbe@0: PG_FREE_IF_COPY(a, 0); \ jbe@0: PG_FREE_IF_COPY(b, 1); \ jbe@0: PG_RETURN_BOOL(res); \ jbe@0: } jbe@0: jbe@0: /* macro for a B-tree comparison function (with detoasting) */ jbe@0: #define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \ jbe@0: PG_FUNCTION_INFO_V1(func); \ jbe@0: Datum func(PG_FUNCTION_ARGS) { \ jbe@0: int32_t res; \ jbe@0: type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ jbe@0: type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ jbe@0: res = cmpfunc(a, b); \ jbe@0: PG_FREE_IF_COPY(a, 0); \ jbe@0: PG_FREE_IF_COPY(b, 1); \ jbe@0: PG_RETURN_INT32(res); \ jbe@0: } jbe@0: jbe@0: /* B-tree operators and comparison function for point */ jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <) jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=) jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==) jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=) jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=) jbe@0: PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >) jbe@0: PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp) jbe@0: jbe@0: /* B-tree operators and comparison function for box */ jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <) jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=) jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==) jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=) jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=) jbe@0: PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >) jbe@0: PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp) jbe@0: jbe@0: /* B-tree operators and comparison function for circle */ jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <) jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=) jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==) jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=) jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=) jbe@0: PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >) jbe@0: PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp) jbe@0: jbe@0: jbe@0: /*--------------------------------* jbe@0: * GiST index support functions * jbe@0: *--------------------------------*/ jbe@0: jbe@0: /* GiST "consistent" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_consistent); jbe@0: Datum pgl_gist_consistent(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); jbe@0: pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); jbe@0: StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); jbe@0: bool *recheck = (bool *)PG_GETARG_POINTER(4); jbe@0: /* demand recheck because index and query methods are lossy */ jbe@0: *recheck = true; jbe@10: /* strategy number aliases for different operators using the same strategy */ jbe@10: strategy %= 100; jbe@0: /* strategy number 11: equality of two points */ jbe@0: if (strategy == 11) { jbe@0: /* query datum is another point */ jbe@0: pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); jbe@0: /* convert other point to key */ jbe@0: pgl_pointkey querykey; jbe@0: pgl_point_to_key(query, querykey); jbe@0: /* return true if both keys overlap */ jbe@0: PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); jbe@0: } jbe@0: /* strategy number 13: equality of two circles */ jbe@0: if (strategy == 13) { jbe@0: /* query datum is another circle */ jbe@0: pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: /* convert other circle to key */ jbe@0: pgl_areakey querykey; jbe@0: pgl_circle_to_key(query, querykey); jbe@0: /* return true if both keys overlap */ jbe@0: PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); jbe@0: } jbe@0: /* for all remaining strategies, keys on empty objects produce no match */ jbe@0: /* (check necessary because query radius may be infinite) */ jbe@0: if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false); jbe@0: /* strategy number 21: overlapping with point */ jbe@0: if (strategy == 21) { jbe@0: /* query datum is a point */ jbe@0: pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); jbe@0: /* return true if estimated distance (allowed to be smaller than real jbe@0: distance) between index key and point is zero */ jbe@0: PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0); jbe@0: } jbe@0: /* strategy number 22: (point) overlapping with box */ jbe@0: if (strategy == 22) { jbe@0: /* query datum is a box */ jbe@0: pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1); jbe@0: /* determine bounding box of indexed key */ jbe@0: pgl_box keybox; jbe@0: pgl_key_to_box(key, &keybox); jbe@0: /* return true if query box overlaps with bounding box of indexed key */ jbe@0: PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox)); jbe@0: } jbe@0: /* strategy number 23: overlapping with circle */ jbe@0: if (strategy == 23) { jbe@0: /* query datum is a circle */ jbe@0: pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: /* return true if estimated distance (allowed to be smaller than real jbe@0: distance) between index key and circle center is smaller than radius */ jbe@0: PG_RETURN_BOOL( jbe@0: pgl_estimate_key_distance(key, &(query->center)) <= query->radius jbe@0: ); jbe@0: } jbe@0: /* strategy number 24: overlapping with cluster */ jbe@0: if (strategy == 24) { jbe@0: bool retval; /* return value */ jbe@0: /* query datum is a cluster */ jbe@0: pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@0: /* return true if estimated distance (allowed to be smaller than real jbe@0: distance) between index key and circle center is smaller than radius */ jbe@0: retval = ( jbe@0: pgl_estimate_key_distance(key, &(query->bounding.center)) <= jbe@0: query->bounding.radius jbe@0: ); jbe@0: PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */ jbe@0: PG_RETURN_BOOL(retval); jbe@0: } jbe@0: /* throw error for any unknown strategy number */ jbe@0: elog(ERROR, "unrecognized strategy number: %d", strategy); jbe@0: } jbe@0: jbe@0: /* GiST "union" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_union); jbe@0: Datum pgl_gist_union(PG_FUNCTION_ARGS) { jbe@0: GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); jbe@0: pgl_keyptr out; /* return value (to be palloc'ed) */ jbe@0: int i; jbe@0: /* determine key size */ jbe@0: size_t keysize = PGL_KEY_IS_AREAKEY( jbe@0: (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key) jbe@0: ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey); jbe@0: /* begin with first key as result */ jbe@0: out = palloc(keysize); jbe@0: memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize); jbe@0: /* unite current result with second, third, etc. key */ jbe@0: for (i=1; in; i++) { jbe@0: pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)); jbe@0: } jbe@0: /* return result */ jbe@0: PG_RETURN_POINTER(out); jbe@0: } jbe@0: jbe@0: /* GiST "compress" support function for indicis on points */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint); jbe@0: Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); jbe@0: GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ jbe@0: /* only transform new leaves */ jbe@0: if (entry->leafkey) { jbe@0: /* get point to be transformed */ jbe@0: pgl_point *point = (pgl_point *)DatumGetPointer(entry->key); jbe@0: /* allocate memory for key */ jbe@0: pgl_keyptr key = palloc(sizeof(pgl_pointkey)); jbe@0: /* transform point to key */ jbe@0: pgl_point_to_key(point, key); jbe@0: /* create new GISTENTRY structure as return value */ jbe@0: retval = palloc(sizeof(GISTENTRY)); jbe@0: gistentryinit( jbe@0: *retval, PointerGetDatum(key), jbe@0: entry->rel, entry->page, entry->offset, FALSE jbe@0: ); jbe@0: } else { jbe@0: /* inner nodes have already been transformed */ jbe@0: retval = entry; jbe@0: } jbe@0: /* return pointer to old or new GISTENTRY structure */ jbe@0: PG_RETURN_POINTER(retval); jbe@0: } jbe@0: jbe@0: /* GiST "compress" support function for indicis on circles */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle); jbe@0: Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); jbe@0: GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ jbe@0: /* only transform new leaves */ jbe@0: if (entry->leafkey) { jbe@0: /* get circle to be transformed */ jbe@0: pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key); jbe@0: /* allocate memory for key */ jbe@0: pgl_keyptr key = palloc(sizeof(pgl_areakey)); jbe@0: /* transform circle to key */ jbe@0: pgl_circle_to_key(circle, key); jbe@0: /* create new GISTENTRY structure as return value */ jbe@0: retval = palloc(sizeof(GISTENTRY)); jbe@0: gistentryinit( jbe@0: *retval, PointerGetDatum(key), jbe@0: entry->rel, entry->page, entry->offset, FALSE jbe@0: ); jbe@0: } else { jbe@0: /* inner nodes have already been transformed */ jbe@0: retval = entry; jbe@0: } jbe@0: /* return pointer to old or new GISTENTRY structure */ jbe@0: PG_RETURN_POINTER(retval); jbe@0: } jbe@0: jbe@0: /* GiST "compress" support function for indices on clusters */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster); jbe@0: Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); jbe@0: GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ jbe@0: /* only transform new leaves */ jbe@0: if (entry->leafkey) { jbe@0: /* get cluster to be transformed (detoasting necessary!) */ jbe@0: pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key); jbe@0: /* allocate memory for key */ jbe@0: pgl_keyptr key = palloc(sizeof(pgl_areakey)); jbe@0: /* transform cluster to key */ jbe@0: pgl_circle_to_key(&(cluster->bounding), key); jbe@0: /* create new GISTENTRY structure as return value */ jbe@0: retval = palloc(sizeof(GISTENTRY)); jbe@0: gistentryinit( jbe@0: *retval, PointerGetDatum(key), jbe@0: entry->rel, entry->page, entry->offset, FALSE jbe@0: ); jbe@0: /* free detoasted datum */ jbe@0: if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster); jbe@0: } else { jbe@0: /* inner nodes have already been transformed */ jbe@0: retval = entry; jbe@0: } jbe@0: /* return pointer to old or new GISTENTRY structure */ jbe@0: PG_RETURN_POINTER(retval); jbe@0: } jbe@0: jbe@0: /* GiST "decompress" support function for indices */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_decompress); jbe@0: Datum pgl_gist_decompress(PG_FUNCTION_ARGS) { jbe@0: /* return passed pointer without transformation */ jbe@0: PG_RETURN_POINTER(PG_GETARG_POINTER(0)); jbe@0: } jbe@0: jbe@0: /* GiST "penalty" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_penalty); jbe@0: Datum pgl_gist_penalty(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0); jbe@0: GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1); jbe@0: float *penalty = (float *)PG_GETARG_POINTER(2); jbe@0: /* get original key and key to insert */ jbe@0: pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key); jbe@0: pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key); jbe@0: /* copy original key */ jbe@0: union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key; jbe@0: if (PGL_KEY_IS_AREAKEY(orig)) { jbe@0: memcpy(union_key.areakey, orig, sizeof(union_key.areakey)); jbe@0: } else { jbe@0: memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey)); jbe@0: } jbe@0: /* calculate union of both keys */ jbe@0: pgl_unite_keys((pgl_keyptr)&union_key, new); jbe@0: /* penalty equal to reduction of key length (logarithm of added area) */ jbe@0: /* (return value by setting referenced value and returning pointer) */ jbe@0: *penalty = ( jbe@0: PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key) jbe@0: ); jbe@0: PG_RETURN_POINTER(penalty); jbe@0: } jbe@0: jbe@0: /* GiST "picksplit" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_picksplit); jbe@0: Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) { jbe@0: GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); jbe@0: GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1); jbe@0: OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */ jbe@0: union { jbe@0: pgl_pointkey pointkey; jbe@0: pgl_areakey areakey; jbe@0: } union_all; /* union of all keys (to be calculated from scratch) jbe@0: (later cut in half) */ jbe@0: int is_areakey = PGL_KEY_IS_AREAKEY( jbe@0: (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key) jbe@0: ); jbe@0: int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey); jbe@0: pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */ jbe@0: pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */ jbe@0: pgl_keyptr key; /* current key to be processed */ jbe@0: /* allocate memory for array of left and right keys, set counts to zero */ jbe@0: v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); jbe@0: v->spl_nleft = 0; jbe@0: v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); jbe@0: v->spl_nright = 0; jbe@0: /* calculate union of all keys from scratch */ jbe@0: memcpy( jbe@0: (pgl_keyptr)&union_all, jbe@0: (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key), jbe@0: keysize jbe@0: ); jbe@0: for (i=FirstOffsetNumber+1; in; i=OffsetNumberNext(i)) { jbe@0: pgl_unite_keys( jbe@0: (pgl_keyptr)&union_all, jbe@0: (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key) jbe@0: ); jbe@0: } jbe@0: /* check if trivial split is necessary due to exhausted key length */ jbe@0: /* (Note: keys for empty objects must have node depth set to maximum) */ jbe@0: if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == ( jbe@0: is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH jbe@0: )) { jbe@0: /* half of all keys go left */ jbe@0: for ( jbe@0: i=FirstOffsetNumber; jbe@0: in - FirstOffsetNumber)/2; jbe@0: i=OffsetNumberNext(i) jbe@0: ) { jbe@0: /* pointer to current key */ jbe@0: key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); jbe@0: /* update unionL */ jbe@0: /* check if key is first key that goes left */ jbe@0: if (!v->spl_nleft) { jbe@0: /* first key that goes left is just copied to unionL */ jbe@0: memcpy(unionL, key, keysize); jbe@0: } else { jbe@0: /* unite current value and next key */ jbe@0: pgl_unite_keys(unionL, key); jbe@0: } jbe@0: /* append offset number to list of keys that go left */ jbe@0: v->spl_left[v->spl_nleft++] = i; jbe@0: } jbe@0: /* other half goes right */ jbe@0: for ( jbe@0: i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; jbe@0: in; jbe@0: i=OffsetNumberNext(i) jbe@0: ) { jbe@0: /* pointer to current key */ jbe@0: key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); jbe@0: /* update unionR */ jbe@0: /* check if key is first key that goes right */ jbe@0: if (!v->spl_nright) { jbe@0: /* first key that goes right is just copied to unionR */ jbe@0: memcpy(unionR, key, keysize); jbe@0: } else { jbe@0: /* unite current value and next key */ jbe@0: pgl_unite_keys(unionR, key); jbe@0: } jbe@0: /* append offset number to list of keys that go right */ jbe@0: v->spl_right[v->spl_nright++] = i; jbe@0: } jbe@0: } jbe@0: /* otherwise, a non-trivial split is possible */ jbe@0: else { jbe@0: /* cut covered area in half */ jbe@0: /* (union_all then refers to area of keys that go left) */ jbe@0: /* check if union of all keys covers empty and non-empty objects */ jbe@0: if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) { jbe@0: /* if yes, split into empty and non-empty objects */ jbe@0: pgl_key_set_empty((pgl_keyptr)&union_all); jbe@0: } else { jbe@0: /* otherwise split by next bit */ jbe@0: ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++; jbe@0: /* NOTE: type bit conserved */ jbe@0: } jbe@0: /* determine for each key if it goes left or right */ jbe@0: for (i=FirstOffsetNumber; in; i=OffsetNumberNext(i)) { jbe@0: /* pointer to current key */ jbe@0: key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); jbe@0: /* keys within one half of the area go left */ jbe@0: if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) { jbe@0: /* update unionL */ jbe@0: /* check if key is first key that goes left */ jbe@0: if (!v->spl_nleft) { jbe@0: /* first key that goes left is just copied to unionL */ jbe@0: memcpy(unionL, key, keysize); jbe@0: } else { jbe@0: /* unite current value of unionL and processed key */ jbe@0: pgl_unite_keys(unionL, key); jbe@0: } jbe@0: /* append offset number to list of keys that go left */ jbe@0: v->spl_left[v->spl_nleft++] = i; jbe@0: } jbe@0: /* the other keys go right */ jbe@0: else { jbe@0: /* update unionR */ jbe@0: /* check if key is first key that goes right */ jbe@0: if (!v->spl_nright) { jbe@0: /* first key that goes right is just copied to unionR */ jbe@0: memcpy(unionR, key, keysize); jbe@0: } else { jbe@0: /* unite current value of unionR and processed key */ jbe@0: pgl_unite_keys(unionR, key); jbe@0: } jbe@0: /* append offset number to list of keys that go right */ jbe@0: v->spl_right[v->spl_nright++] = i; jbe@0: } jbe@0: } jbe@0: } jbe@0: /* store unions in return value */ jbe@0: v->spl_ldatum = PointerGetDatum(unionL); jbe@0: v->spl_rdatum = PointerGetDatum(unionR); jbe@0: /* return all results */ jbe@0: PG_RETURN_POINTER(v); jbe@0: } jbe@0: jbe@0: /* GiST "same"/"equal" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_same); jbe@0: Datum pgl_gist_same(PG_FUNCTION_ARGS) { jbe@0: pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0); jbe@0: pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1); jbe@0: bool *boolptr = (bool *)PG_GETARG_POINTER(2); jbe@0: /* two keys are equal if they are binary equal */ jbe@0: /* (return result by setting referenced boolean and returning pointer) */ jbe@0: *boolptr = !memcmp( jbe@0: key1, jbe@0: key2, jbe@0: PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey) jbe@0: ); jbe@0: PG_RETURN_POINTER(boolptr); jbe@0: } jbe@0: jbe@0: /* GiST "distance" support function */ jbe@0: PG_FUNCTION_INFO_V1(pgl_gist_distance); jbe@0: Datum pgl_gist_distance(PG_FUNCTION_ARGS) { jbe@0: GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0); jbe@0: pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); jbe@0: StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); jbe@0: bool *recheck = (bool *)PG_GETARG_POINTER(4); jbe@0: double distance; /* return value */ jbe@0: /* demand recheck because distance is just an estimation */ jbe@0: /* (real distance may be bigger) */ jbe@0: *recheck = true; jbe@10: /* strategy number aliases for different operators using the same strategy */ jbe@10: strategy %= 100; jbe@0: /* strategy number 31: distance to point */ jbe@0: if (strategy == 31) { jbe@0: /* query datum is a point */ jbe@0: pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); jbe@0: /* use pgl_estimate_pointkey_distance() function to compute result */ jbe@0: distance = pgl_estimate_key_distance(key, query); jbe@0: /* avoid infinity (reserved!) */ jbe@0: if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; jbe@0: /* return result */ jbe@0: PG_RETURN_FLOAT8(distance); jbe@0: } jbe@0: /* strategy number 33: distance to circle */ jbe@0: if (strategy == 33) { jbe@0: /* query datum is a circle */ jbe@0: pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); jbe@0: /* estimate distance to circle center and substract circle radius */ jbe@0: distance = ( jbe@0: pgl_estimate_key_distance(key, &(query->center)) - query->radius jbe@0: ); jbe@0: /* convert non-positive values to zero and avoid infinity (reserved!) */ jbe@0: if (distance <= 0) distance = 0; jbe@0: else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; jbe@0: /* return result */ jbe@0: PG_RETURN_FLOAT8(distance); jbe@0: } jbe@0: /* strategy number 34: distance to cluster */ jbe@0: if (strategy == 34) { jbe@0: /* query datum is a cluster */ jbe@0: pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); jbe@0: /* estimate distance to bounding center and substract bounding radius */ jbe@0: distance = ( jbe@0: pgl_estimate_key_distance(key, &(query->bounding.center)) - jbe@0: query->bounding.radius jbe@0: ); jbe@0: /* convert non-positive values to zero and avoid infinity (reserved!) */ jbe@0: if (distance <= 0) distance = 0; jbe@0: else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; jbe@0: /* free detoasted cluster (if copy) */ jbe@0: PG_FREE_IF_COPY(query, 1); jbe@0: /* return result */ jbe@0: PG_RETURN_FLOAT8(distance); jbe@0: } jbe@0: /* throw error for any unknown strategy number */ jbe@0: elog(ERROR, "unrecognized strategy number: %d", strategy); jbe@0: } jbe@0: