pgLatLon is a spatial database extension for the PostgreSQL object-relational database management system providing geographic data types and spatial indexing @@ -40,8 +40,8 @@
It is also possible to compile and install the extension without GNU Make as follows:
-cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c
-cp latlon-v0007.so `pg_config --pkglibdir`
+cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c
+cp latlon-v0008.so `pg_config --pkglibdir`
cp latlon.control `pg_config --sharedir`/extension/
cp latlon--*.sql `pg_config --sharedir`/extension/
@@ -486,6 +486,40 @@
Same as epoint(float8, float8) but with arguments reversed.
+fair_distance(ecluster, epoint, samples int4 = 10000)
+
+When working with user-generated content, users may be tempted to create
+intentionally oversized objects in order to optimize search results in an
+unfair manner. The fair_distance function aims to handle this by returning an
+adjusted distance (i.e. distance increased by a penalty) if a geographic object
+(the ecluster) consists of more than one point.
+
+The first argument to this function is an ecluster, the second argument is a
+search point (epoint), and the third argument is an interger related to the
+precision (higher precision will require more computation time).
+
+The penalty by which the returned distance is increased fulfills (at least) the
+following properties:
+
+
+- For search points far away from the
ecluster (i.e. distance approaching
+infinity), the penalty approaches zero (i.e. fair_distance behaves the same
+as distance).
+- If the
ecluster consists of a set of points, the penalty for a search point
+close to one of that points (closer than half of the minimum distance between
+each pair of points in the ecluster) is chosen in such a way that the
+adjusted distance is equal to the distance from the search point to the
+closest point in the ecluster multiplied by the square root of the count of
+points in the ecluster.
+
+
+The function interally uses a Monte Carlo simulation to compute the result. The
+third parameter (which defaults to 10000) can be used to adjust the number of
+samples taken. It is ensured that the penalty is always positive, i.e. results
+returned by the fair_distance function are always equal to or greater than
+the results returned by the distance function regardless of stochastic
+effects.
+
GeoJSON_to_epoint(jsonb, text)
Maps a GeoJSON object of type "Point" or "Feature" (which contains a
diff -r a90c0d750405 -r 1b9cd45e9e48 README.mkd
--- a/README.mkd Fri Oct 21 12:57:46 2016 +0200
+++ b/README.mkd Tue Oct 25 18:44:43 2016 +0200
@@ -1,4 +1,4 @@
-pgLatLon v0.8 documentation
+pgLatLon v0.9 documentation
===========================
pgLatLon is a spatial database extension for the PostgreSQL object-relational
@@ -39,8 +39,8 @@
It is also possible to compile and install the extension without GNU Make as
follows:
- cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c
- cp latlon-v0007.so `pg_config --pkglibdir`
+ cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c
+ cp latlon-v0008.so `pg_config --pkglibdir`
cp latlon.control `pg_config --sharedir`/extension/
cp latlon--*.sql `pg_config --sharedir`/extension/
@@ -469,6 +469,38 @@
Same as `epoint(float8, float8)` but with arguments reversed.
+#### `fair_distance(ecluster, epoint,` samples `int4 = 10000)`
+
+When working with user-generated content, users may be tempted to create
+intentionally oversized objects in order to optimize search results in an
+unfair manner. The `fair_distance` function aims to handle this by returning an
+adjusted distance (i.e. distance increased by a penalty) if a geographic object
+(the `ecluster`) consists of more than one point.
+
+The first argument to this function is an `ecluster`, the second argument is a
+search point (`epoint`), and the third argument is an interger related to the
+precision (higher precision will require more computation time).
+
+The penalty by which the returned distance is increased fulfills (at least) the
+following properties:
+
+* For search points far away from the `ecluster` (i.e. distance approaching
+ infinity), the penalty approaches zero (i.e. `fair_distance` behaves the same
+ as `distance`).
+* If the `ecluster` consists of a set of points, the penalty for a search point
+ close to one of that points (closer than half of the minimum distance between
+ each pair of points in the `ecluster`) is chosen in such a way that the
+ adjusted distance is equal to the distance from the search point to the
+ closest point in the `ecluster` multiplied by the square root of the count of
+ points in the `ecluster`.
+
+The function interally uses a Monte Carlo simulation to compute the result. The
+third parameter (which defaults to 10000) can be used to adjust the number of
+samples taken. It is ensured that the penalty is always positive, i.e. results
+returned by the `fair_distance` function are always equal to or greater than
+the results returned by the `distance` function regardless of stochastic
+effects.
+
#### `GeoJSON_to_epoint(jsonb, text)`
Maps a GeoJSON object of type "Point" or "Feature" (which contains a
diff -r a90c0d750405 -r 1b9cd45e9e48 latlon--0.9.sql
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/latlon--0.9.sql Tue Oct 25 18:44:43 2016 +0200
@@ -0,0 +1,1692 @@
+
+----------------------------------------
+-- forward declarations (shell types) --
+----------------------------------------
+
+CREATE TYPE epoint;
+CREATE TYPE ebox;
+CREATE TYPE ecircle;
+CREATE TYPE ecluster;
+
+
+------------------------------------------------------------
+-- dummy input/output functions for dummy index key types --
+------------------------------------------------------------
+
+CREATE FUNCTION ekey_point_in_dummy(cstring)
+ RETURNS ekey_point
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_notimpl';
+
+CREATE FUNCTION ekey_point_out_dummy(ekey_point)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_notimpl';
+
+CREATE FUNCTION ekey_area_in_dummy(cstring)
+ RETURNS ekey_area
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_notimpl';
+
+CREATE FUNCTION ekey_area_out_dummy(ekey_area)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_notimpl';
+
+
+--------------------------
+-- text input functions --
+--------------------------
+
+CREATE FUNCTION epoint_in(cstring)
+ RETURNS epoint
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_in';
+
+CREATE FUNCTION ebox_in(cstring)
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_in';
+
+CREATE FUNCTION ecircle_in(cstring)
+ RETURNS ecircle
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_in';
+
+CREATE FUNCTION ecluster_in(cstring)
+ RETURNS ecluster
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_in';
+
+
+---------------------------
+-- text output functions --
+---------------------------
+
+CREATE FUNCTION epoint_out(epoint)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_out';
+
+CREATE FUNCTION ebox_out(ebox)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_out';
+
+CREATE FUNCTION ecircle_out(ecircle)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_out';
+
+CREATE FUNCTION ecluster_out(ecluster)
+ RETURNS cstring
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_out';
+
+
+--------------------------
+-- binary I/O functions --
+--------------------------
+
+CREATE FUNCTION epoint_recv(internal)
+ RETURNS epoint
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_recv';
+
+CREATE FUNCTION ebox_recv(internal)
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_recv';
+
+CREATE FUNCTION ecircle_recv(internal)
+ RETURNS ecircle
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_recv';
+
+CREATE FUNCTION epoint_send(epoint)
+ RETURNS bytea
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_send';
+
+CREATE FUNCTION ebox_send(ebox)
+ RETURNS bytea
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_send';
+
+CREATE FUNCTION ecircle_send(ecircle)
+ RETURNS bytea
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_send';
+
+
+-----------------------------------------------
+-- type definitions of dummy index key types --
+-----------------------------------------------
+
+CREATE TYPE ekey_point (
+ internallength = 8,
+ input = ekey_point_in_dummy,
+ output = ekey_point_out_dummy,
+ alignment = char );
+
+CREATE TYPE ekey_area (
+ internallength = 9,
+ input = ekey_area_in_dummy,
+ output = ekey_area_out_dummy,
+ alignment = char );
+
+
+------------------------------------------
+-- definitions of geographic data types --
+------------------------------------------
+
+CREATE TYPE epoint (
+ internallength = 16,
+ input = epoint_in,
+ output = epoint_out,
+ receive = epoint_recv,
+ send = epoint_send,
+ alignment = double );
+
+CREATE TYPE ebox (
+ internallength = 32,
+ input = ebox_in,
+ output = ebox_out,
+ receive = ebox_recv,
+ send = ebox_send,
+ alignment = double );
+
+CREATE TYPE ecircle (
+ internallength = 24,
+ input = ecircle_in,
+ output = ecircle_out,
+ receive = ecircle_recv,
+ send = ecircle_send,
+ alignment = double );
+
+CREATE TYPE ecluster (
+ internallength = VARIABLE,
+ input = ecluster_in,
+ output = ecluster_out,
+ alignment = double,
+ storage = external );
+
+
+--------------------
+-- B-tree support --
+--------------------
+
+-- begin of B-tree support for epoint
+
+CREATE FUNCTION epoint_btree_lt(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_lt';
+
+CREATE FUNCTION epoint_btree_le(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_le';
+
+CREATE FUNCTION epoint_btree_eq(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_eq';
+
+CREATE FUNCTION epoint_btree_ne(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ne';
+
+CREATE FUNCTION epoint_btree_ge(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ge';
+
+CREATE FUNCTION epoint_btree_gt(epoint, epoint)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_gt';
+
+CREATE OPERATOR <<< (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_lt,
+ commutator = >>>,
+ negator = >>>=,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR <<<= (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_le,
+ commutator = >>>=,
+ negator = >>>,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR = (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_eq,
+ commutator = =,
+ negator = <>,
+ restrict = eqsel,
+ join = eqjoinsel,
+ merges
+);
+
+CREATE OPERATOR <> (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_eq,
+ commutator = <>,
+ negator = =,
+ restrict = neqsel,
+ join = neqjoinsel
+);
+
+CREATE OPERATOR >>>= (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_ge,
+ commutator = <<<=,
+ negator = <<<,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE OPERATOR >>> (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_btree_gt,
+ commutator = <<<,
+ negator = <<<=,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE FUNCTION epoint_btree_cmp(epoint, epoint)
+ RETURNS int4
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_epoint_cmp';
+
+CREATE OPERATOR CLASS epoint_btree_ops
+ DEFAULT FOR TYPE epoint USING btree AS
+ OPERATOR 1 <<< ,
+ OPERATOR 2 <<<= ,
+ OPERATOR 3 = ,
+ OPERATOR 4 >>>= ,
+ OPERATOR 5 >>> ,
+ FUNCTION 1 epoint_btree_cmp(epoint, epoint);
+
+-- end of B-tree support for epoint
+
+-- begin of B-tree support for ebox
+
+CREATE FUNCTION ebox_btree_lt(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_lt';
+
+CREATE FUNCTION ebox_btree_le(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_le';
+
+CREATE FUNCTION ebox_btree_eq(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_eq';
+
+CREATE FUNCTION ebox_btree_ne(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ne';
+
+CREATE FUNCTION ebox_btree_ge(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ge';
+
+CREATE FUNCTION ebox_btree_gt(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_gt';
+
+CREATE OPERATOR <<< (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_lt,
+ commutator = >>>,
+ negator = >>>=,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR <<<= (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_le,
+ commutator = >>>=,
+ negator = >>>,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR = (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_eq,
+ commutator = =,
+ negator = <>,
+ restrict = eqsel,
+ join = eqjoinsel,
+ merges
+);
+
+CREATE OPERATOR <> (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_eq,
+ commutator = <>,
+ negator = =,
+ restrict = neqsel,
+ join = neqjoinsel
+);
+
+CREATE OPERATOR >>>= (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_ge,
+ commutator = <<<=,
+ negator = <<<,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE OPERATOR >>> (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_btree_gt,
+ commutator = <<<,
+ negator = <<<=,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE FUNCTION ebox_btree_cmp(ebox, ebox)
+ RETURNS int4
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ebox_cmp';
+
+CREATE OPERATOR CLASS ebox_btree_ops
+ DEFAULT FOR TYPE ebox USING btree AS
+ OPERATOR 1 <<< ,
+ OPERATOR 2 <<<= ,
+ OPERATOR 3 = ,
+ OPERATOR 4 >>>= ,
+ OPERATOR 5 >>> ,
+ FUNCTION 1 ebox_btree_cmp(ebox, ebox);
+
+-- end of B-tree support for ebox
+
+-- begin of B-tree support for ecircle
+
+CREATE FUNCTION ecircle_btree_lt(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_lt';
+
+CREATE FUNCTION ecircle_btree_le(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_le';
+
+CREATE FUNCTION ecircle_btree_eq(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_eq';
+
+CREATE FUNCTION ecircle_btree_ne(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ne';
+
+CREATE FUNCTION ecircle_btree_ge(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ge';
+
+CREATE FUNCTION ecircle_btree_gt(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_gt';
+
+CREATE OPERATOR <<< (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_lt,
+ commutator = >>>,
+ negator = >>>=,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR <<<= (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_le,
+ commutator = >>>=,
+ negator = >>>,
+ restrict = scalarltsel,
+ join = scalarltjoinsel
+);
+
+CREATE OPERATOR = (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_eq,
+ commutator = =,
+ negator = <>,
+ restrict = eqsel,
+ join = eqjoinsel,
+ merges
+);
+
+CREATE OPERATOR <> (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_eq,
+ commutator = <>,
+ negator = =,
+ restrict = neqsel,
+ join = neqjoinsel
+);
+
+CREATE OPERATOR >>>= (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_ge,
+ commutator = <<<=,
+ negator = <<<,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE OPERATOR >>> (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_btree_gt,
+ commutator = <<<,
+ negator = <<<=,
+ restrict = scalargtsel,
+ join = scalargtjoinsel
+);
+
+CREATE FUNCTION ecircle_btree_cmp(ecircle, ecircle)
+ RETURNS int4
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_cmp';
+
+CREATE OPERATOR CLASS ecircle_btree_ops
+ DEFAULT FOR TYPE ecircle USING btree AS
+ OPERATOR 1 <<< ,
+ OPERATOR 2 <<<= ,
+ OPERATOR 3 = ,
+ OPERATOR 4 >>>= ,
+ OPERATOR 5 >>> ,
+ FUNCTION 1 ecircle_btree_cmp(ecircle, ecircle);
+
+-- end of B-tree support for ecircle
+
+
+----------------
+-- type casts --
+----------------
+
+CREATE FUNCTION cast_epoint_to_ebox(epoint)
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_to_ebox';
+
+CREATE CAST (epoint AS ebox) WITH FUNCTION cast_epoint_to_ebox(epoint);
+
+CREATE FUNCTION cast_epoint_to_ecircle(epoint)
+ RETURNS ecircle
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecircle';
+
+CREATE CAST (epoint AS ecircle) WITH FUNCTION cast_epoint_to_ecircle(epoint);
+
+CREATE FUNCTION cast_epoint_to_ecluster(epoint)
+ RETURNS ecluster
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecluster';
+
+CREATE CAST (epoint AS ecluster) WITH FUNCTION cast_epoint_to_ecluster(epoint);
+
+CREATE FUNCTION cast_ebox_to_ecluster(ebox)
+ RETURNS ecluster
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_to_ecluster';
+
+CREATE CAST (ebox AS ecluster) WITH FUNCTION cast_ebox_to_ecluster(ebox);
+
+
+---------------------------
+-- constructor functions --
+---------------------------
+
+CREATE FUNCTION epoint(float8, float8)
+ RETURNS epoint
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_epoint';
+
+CREATE FUNCTION epoint_latlon(float8, float8)
+ RETURNS epoint
+ LANGUAGE SQL IMMUTABLE STRICT AS $$
+ SELECT epoint($1, $2)
+ $$;
+
+CREATE FUNCTION epoint_lonlat(float8, float8)
+ RETURNS epoint
+ LANGUAGE SQL IMMUTABLE STRICT AS $$
+ SELECT epoint($2, $1)
+ $$;
+
+CREATE FUNCTION empty_ebox()
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_empty_ebox';
+
+CREATE FUNCTION ebox(float8, float8, float8, float8)
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_ebox';
+
+CREATE FUNCTION ebox(epoint, epoint)
+ RETURNS ebox
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_ebox_from_epoints';
+
+CREATE FUNCTION ecircle(float8, float8, float8)
+ RETURNS ecircle
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_ecircle';
+
+CREATE FUNCTION ecircle(epoint, float8)
+ RETURNS ecircle
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_create_ecircle_from_epoint';
+
+CREATE FUNCTION ecluster_concat(ecluster[])
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT array_to_string($1, ' ')::ecluster
+ $$;
+
+CREATE FUNCTION ecluster_concat(ecluster, ecluster)
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT ($1::text || ' ' || $2::text)::ecluster
+ $$;
+
+CREATE FUNCTION ecluster_create_multipoint(epoint[])
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT
+ array_to_string(array_agg('point (' || unnest || ')'), ' ')::ecluster
+ FROM unnest($1)
+ $$;
+
+CREATE FUNCTION ecluster_create_path(epoint[])
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
+ ('path (' || array_to_string($1, ' ') || ')')::ecluster
+ END
+ FROM array_to_string($1, ' ') AS "str"
+ $$;
+
+CREATE FUNCTION ecluster_create_outline(epoint[])
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
+ ('outline (' || array_to_string($1, ' ') || ')')::ecluster
+ END
+ FROM array_to_string($1, ' ') AS "str"
+ $$;
+
+CREATE FUNCTION ecluster_create_polygon(epoint[])
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
+ ('polygon (' || array_to_string($1, ' ') || ')')::ecluster
+ END
+ FROM array_to_string($1, ' ') AS "str"
+ $$;
+
+
+----------------------
+-- getter functions --
+----------------------
+
+CREATE FUNCTION latitude(epoint)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_lat';
+
+CREATE FUNCTION longitude(epoint)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_lon';
+
+CREATE FUNCTION min_latitude(ebox)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_lat_min';
+
+CREATE FUNCTION max_latitude(ebox)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_lat_max';
+
+CREATE FUNCTION min_longitude(ebox)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_lon_min';
+
+CREATE FUNCTION max_longitude(ebox)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_lon_max';
+
+CREATE FUNCTION center(ecircle)
+ RETURNS epoint
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_center';
+
+CREATE FUNCTION radius(ecircle)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_radius';
+
+CREATE FUNCTION ecluster_extract_points(ecluster)
+ RETURNS SETOF epoint
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT "match"[2]::epoint
+ FROM regexp_matches($1::text, e'(^| )point \\(([^)]+)\\)', 'g') AS "match"
+ $$;
+
+CREATE FUNCTION ecluster_extract_paths(ecluster)
+ RETURNS SETOF epoint[]
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT (
+ SELECT array_agg("m2"[1]::epoint)
+ FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
+ )
+ FROM regexp_matches($1::text, e'(^| )path \\(([^)]+)\\)', 'g') AS "m1"
+ $$;
+
+CREATE FUNCTION ecluster_extract_outlines(ecluster)
+ RETURNS SETOF epoint[]
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT (
+ SELECT array_agg("m2"[1]::epoint)
+ FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
+ )
+ FROM regexp_matches($1::text, e'(^| )outline \\(([^)]+)\\)', 'g') AS "m1"
+ $$;
+
+CREATE FUNCTION ecluster_extract_polygons(ecluster)
+ RETURNS SETOF epoint[]
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT (
+ SELECT array_agg("m2"[1]::epoint)
+ FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
+ )
+ FROM regexp_matches($1::text, e'(^| )polygon \\(([^)]+)\\)', 'g') AS "m1"
+ $$;
+
+
+---------------
+-- operators --
+---------------
+
+CREATE FUNCTION epoint_ebox_overlap_proc(epoint, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ebox_overlap';
+
+CREATE FUNCTION epoint_ecircle_overlap_proc(epoint, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_overlap';
+
+CREATE FUNCTION epoint_ecluster_overlap_proc(epoint, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_overlap';
+
+CREATE FUNCTION epoint_ecluster_may_overlap_proc(epoint, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_may_overlap';
+
+CREATE FUNCTION ebox_overlap_proc(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_overlap';
+
+CREATE FUNCTION ebox_ecircle_may_overlap_proc(ebox, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_ecircle_may_overlap';
+
+CREATE FUNCTION ebox_ecluster_may_overlap_proc(ebox, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ebox_ecluster_may_overlap';
+
+CREATE FUNCTION ecircle_overlap_proc(ecircle, ecircle)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_overlap';
+
+CREATE FUNCTION ecircle_ecluster_overlap_proc(ecircle, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_overlap';
+
+CREATE FUNCTION ecircle_ecluster_may_overlap_proc(ecircle, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_may_overlap';
+
+CREATE FUNCTION ecluster_overlap_proc(ecluster, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_overlap';
+
+CREATE FUNCTION ecluster_may_overlap_proc(ecluster, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_may_overlap';
+
+CREATE FUNCTION ecluster_contains_proc(ecluster, ecluster)
+ RETURNS boolean
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_contains';
+
+CREATE FUNCTION epoint_distance_proc(epoint, epoint)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_distance';
+
+CREATE FUNCTION epoint_ecircle_distance_proc(epoint, ecircle)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_distance';
+
+CREATE FUNCTION epoint_ecluster_distance_proc(epoint, ecluster)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_distance';
+
+CREATE FUNCTION ecircle_distance_proc(ecircle, ecircle)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_distance';
+
+CREATE FUNCTION ecircle_ecluster_distance_proc(ecircle, ecluster)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_distance';
+
+CREATE FUNCTION ecluster_distance_proc(ecluster, ecluster)
+ RETURNS float8
+ LANGUAGE C IMMUTABLE STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_distance';
+
+CREATE OPERATOR && (
+ leftarg = epoint,
+ rightarg = ebox,
+ procedure = epoint_ebox_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION epoint_ebox_overlap_commutator(ebox, epoint)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
+
+CREATE OPERATOR && (
+ leftarg = ebox,
+ rightarg = epoint,
+ procedure = epoint_ebox_overlap_commutator,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = epoint,
+ rightarg = ecircle,
+ procedure = epoint_ecircle_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION epoint_ecircle_overlap_commutator(ecircle, epoint)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
+
+CREATE OPERATOR && (
+ leftarg = ecircle,
+ rightarg = epoint,
+ procedure = epoint_ecircle_overlap_commutator,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = epoint,
+ rightarg = ecluster,
+ procedure = epoint_ecluster_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION epoint_ecluster_overlap_commutator(ecluster, epoint)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
+
+CREATE OPERATOR && (
+ leftarg = ecluster,
+ rightarg = epoint,
+ procedure = epoint_ecluster_overlap_commutator,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = ecircle,
+ rightarg = ecluster,
+ procedure = ecircle_ecluster_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ecircle_ecluster_overlap_commutator(ecluster, ecircle)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
+
+CREATE OPERATOR && (
+ leftarg = ecluster,
+ rightarg = ecircle,
+ procedure = ecircle_ecluster_overlap_commutator,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR && (
+ leftarg = ecluster,
+ rightarg = ecluster,
+ procedure = ecluster_overlap_proc,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecircle_overlap_castwrap(ebox, ecircle)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
+
+CREATE OPERATOR && (
+ leftarg = ebox,
+ rightarg = ecircle,
+ procedure = ebox_ecircle_overlap_castwrap,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecircle_overlap_castwrap(ecircle, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
+
+CREATE OPERATOR && (
+ leftarg = ecircle,
+ rightarg = ebox,
+ procedure = ebox_ecircle_overlap_castwrap,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecluster_overlap_castwrap(ebox, ecluster)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
+
+CREATE OPERATOR && (
+ leftarg = ebox,
+ rightarg = ecluster,
+ procedure = ebox_ecluster_overlap_castwrap,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecluster_overlap_castwrap(ecluster, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
+
+CREATE OPERATOR && (
+ leftarg = ecluster,
+ rightarg = ebox,
+ procedure = ebox_ecluster_overlap_castwrap,
+ commutator = &&,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR &&+ (
+ leftarg = epoint,
+ rightarg = ecluster,
+ procedure = epoint_ecluster_may_overlap_proc,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION epoint_ecluster_may_overlap_commutator(ecluster, epoint)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
+
+CREATE OPERATOR &&+ (
+ leftarg = ecluster,
+ rightarg = epoint,
+ procedure = epoint_ecluster_may_overlap_commutator,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR &&+ (
+ leftarg = ebox,
+ rightarg = ecircle,
+ procedure = ebox_ecircle_may_overlap_proc,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecircle_may_overlap_commutator(ecircle, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
+
+CREATE OPERATOR &&+ (
+ leftarg = ecircle,
+ rightarg = ebox,
+ procedure = ebox_ecircle_may_overlap_commutator,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR &&+ (
+ leftarg = ebox,
+ rightarg = ecluster,
+ procedure = ebox_ecluster_may_overlap_proc,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecluster_may_overlap_commutator(ecluster, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
+
+CREATE OPERATOR &&+ (
+ leftarg = ecluster,
+ rightarg = ebox,
+ procedure = ebox_ecluster_may_overlap_commutator,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR &&+ (
+ leftarg = ecircle,
+ rightarg = ecluster,
+ procedure = ecircle_ecluster_may_overlap_proc,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ecircle_ecluster_may_overlap_commutator(ecluster, ecircle)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
+
+CREATE OPERATOR &&+ (
+ leftarg = ecluster,
+ rightarg = ecircle,
+ procedure = ecircle_ecluster_may_overlap_commutator,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR &&+ (
+ leftarg = ecluster,
+ rightarg = ecluster,
+ procedure = ecluster_may_overlap_proc,
+ commutator = &&+,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR @> (
+ leftarg = ebox,
+ rightarg = epoint,
+ procedure = epoint_ebox_overlap_commutator,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR <@ (
+ leftarg = epoint,
+ rightarg = ebox,
+ procedure = epoint_ebox_overlap_proc,
+ commutator = @>,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR @> (
+ leftarg = ecluster,
+ rightarg = epoint,
+ procedure = epoint_ecluster_overlap_commutator,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR <@ (
+ leftarg = epoint,
+ rightarg = ecluster,
+ procedure = epoint_ecluster_overlap_proc,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR @> (
+ leftarg = ecluster,
+ rightarg = ecluster,
+ procedure = ecluster_contains_proc,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ecluster_contains_commutator(ecluster, ecluster)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1';
+
+CREATE OPERATOR <@ (
+ leftarg = ecluster,
+ rightarg = ecluster,
+ procedure = ecluster_contains_commutator,
+ commutator = @>,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_contains_castwrap(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2::ecluster';
+
+CREATE OPERATOR @> (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_contains_castwrap,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_contains_swapped_castwrap(ebox, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1::ecluster';
+
+CREATE OPERATOR <@ (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_contains_swapped_castwrap,
+ commutator = @>,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecluster_contains_castwrap(ebox, ecluster)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2';
+
+CREATE OPERATOR @> (
+ leftarg = ebox,
+ rightarg = ecluster,
+ procedure = ebox_ecluster_contains_castwrap,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ebox_ecluster_contains_castwrap(ecluster, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1';
+
+CREATE OPERATOR <@ (
+ leftarg = ecluster,
+ rightarg = ebox,
+ procedure = ebox_ecluster_contains_castwrap,
+ commutator = @>,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ecluster_ebox_contains_castwrap(ecluster, ebox)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 @> $2::ecluster';
+
+CREATE OPERATOR @> (
+ leftarg = ecluster,
+ rightarg = ebox,
+ procedure = ecluster_ebox_contains_castwrap,
+ commutator = <@,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE FUNCTION ecluster_ebox_contains_castwrap(ebox, ecluster)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1::ecluster';
+
+CREATE OPERATOR <@ (
+ leftarg = ebox,
+ rightarg = ecluster,
+ procedure = ecluster_ebox_contains_castwrap,
+ commutator = @>,
+ restrict = areasel,
+ join = areajoinsel
+);
+
+CREATE OPERATOR <-> (
+ leftarg = epoint,
+ rightarg = epoint,
+ procedure = epoint_distance_proc,
+ commutator = <->
+);
+
+CREATE OPERATOR <-> (
+ leftarg = epoint,
+ rightarg = ecircle,
+ procedure = epoint_ecircle_distance_proc,
+ commutator = <->
+);
+
+CREATE FUNCTION epoint_ecircle_distance_commutator(ecircle, epoint)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
+
+CREATE OPERATOR <-> (
+ leftarg = ecircle,
+ rightarg = epoint,
+ procedure = epoint_ecircle_distance_commutator,
+ commutator = <->
+);
+
+CREATE OPERATOR <-> (
+ leftarg = epoint,
+ rightarg = ecluster,
+ procedure = epoint_ecluster_distance_proc,
+ commutator = <->
+);
+
+CREATE FUNCTION epoint_ecluster_distance_commutator(ecluster, epoint)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
+
+CREATE OPERATOR <-> (
+ leftarg = ecluster,
+ rightarg = epoint,
+ procedure = epoint_ecluster_distance_commutator,
+ commutator = <->
+);
+
+CREATE OPERATOR <-> (
+ leftarg = ecircle,
+ rightarg = ecircle,
+ procedure = ecircle_distance_proc,
+ commutator = <->
+);
+
+CREATE OPERATOR <-> (
+ leftarg = ecircle,
+ rightarg = ecluster,
+ procedure = ecircle_ecluster_distance_proc,
+ commutator = <->
+);
+
+CREATE FUNCTION ecircle_ecluster_distance_commutator(ecluster, ecircle)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
+
+CREATE OPERATOR <-> (
+ leftarg = ecluster,
+ rightarg = ecircle,
+ procedure = ecircle_ecluster_distance_commutator,
+ commutator = <->
+);
+
+CREATE OPERATOR <-> (
+ leftarg = ecluster,
+ rightarg = ecluster,
+ procedure = ecluster_distance_proc,
+ commutator = <->
+);
+
+CREATE FUNCTION epoint_ebox_distance_castwrap(epoint, ebox)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
+
+CREATE OPERATOR <-> (
+ leftarg = epoint,
+ rightarg = ebox,
+ procedure = epoint_ebox_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION epoint_ebox_distance_castwrap(ebox, epoint)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
+
+CREATE OPERATOR <-> (
+ leftarg = ebox,
+ rightarg = epoint,
+ procedure = epoint_ebox_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION ebox_distance_castwrap(ebox, ebox)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2::ecluster';
+
+CREATE OPERATOR <-> (
+ leftarg = ebox,
+ rightarg = ebox,
+ procedure = ebox_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION ebox_ecircle_distance_castwrap(ebox, ecircle)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
+
+CREATE OPERATOR <-> (
+ leftarg = ebox,
+ rightarg = ecircle,
+ procedure = ebox_ecircle_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION ebox_ecircle_distance_castwrap(ecircle, ebox)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
+
+CREATE OPERATOR <-> (
+ leftarg = ecircle,
+ rightarg = ebox,
+ procedure = ebox_ecircle_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION ebox_ecluster_distance_castwrap(ebox, ecluster)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
+
+CREATE OPERATOR <-> (
+ leftarg = ebox,
+ rightarg = ecluster,
+ procedure = ebox_ecluster_distance_castwrap,
+ commutator = <->
+);
+
+CREATE FUNCTION ebox_ecluster_distance_castwrap(ecluster, ebox)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
+
+CREATE OPERATOR <-> (
+ leftarg = ecluster,
+ rightarg = ebox,
+ procedure = ebox_ecluster_distance_castwrap,
+ commutator = <->
+);
+
+
+---------------------
+-- other functions --
+---------------------
+
+CREATE FUNCTION monte_carlo_area(ecluster, int4 = 10000)
+ RETURNS float8
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_monte_carlo_area';
+
+CREATE FUNCTION fair_distance(ecluster, epoint, int4 = 10000)
+ RETURNS float8
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_ecluster_epoint_fair_distance';
+
+
+----------------
+-- GiST index --
+----------------
+
+CREATE FUNCTION pgl_gist_consistent(internal, internal, smallint, oid, internal)
+ RETURNS boolean
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_consistent';
+
+CREATE FUNCTION pgl_gist_union(internal, internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_union';
+
+CREATE FUNCTION pgl_gist_compress_epoint(internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_compress_epoint';
+
+CREATE FUNCTION pgl_gist_compress_ecircle(internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecircle';
+
+CREATE FUNCTION pgl_gist_compress_ecluster(internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecluster';
+
+CREATE FUNCTION pgl_gist_decompress(internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_decompress';
+
+CREATE FUNCTION pgl_gist_penalty(internal, internal, internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_penalty';
+
+CREATE FUNCTION pgl_gist_picksplit(internal, internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_picksplit';
+
+CREATE FUNCTION pgl_gist_same(internal, internal, internal)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_same';
+
+CREATE FUNCTION pgl_gist_distance(internal, internal, smallint, oid)
+ RETURNS internal
+ LANGUAGE C STRICT
+ AS '$libdir/latlon-v0008', 'pgl_gist_distance';
+
+CREATE OPERATOR CLASS epoint_ops
+ DEFAULT FOR TYPE epoint USING gist AS
+ OPERATOR 11 = ,
+ OPERATOR 22 && (epoint, ebox),
+ OPERATOR 222 <@ (epoint, ebox),
+ OPERATOR 23 && (epoint, ecircle),
+ OPERATOR 24 && (epoint, ecluster),
+ OPERATOR 124 &&+ (epoint, ecluster),
+ OPERATOR 224 <@ (epoint, ecluster),
+ OPERATOR 31 <-> (epoint, epoint) FOR ORDER BY float_ops,
+ OPERATOR 32 <-> (epoint, ebox) FOR ORDER BY float_ops,
+ OPERATOR 33 <-> (epoint, ecircle) FOR ORDER BY float_ops,
+ OPERATOR 34 <-> (epoint, ecluster) FOR ORDER BY float_ops,
+ FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
+ FUNCTION 2 pgl_gist_union(internal, internal),
+ FUNCTION 3 pgl_gist_compress_epoint(internal),
+ FUNCTION 4 pgl_gist_decompress(internal),
+ FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
+ FUNCTION 6 pgl_gist_picksplit(internal, internal),
+ FUNCTION 7 pgl_gist_same(internal, internal, internal),
+ FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
+ STORAGE ekey_point;
+
+CREATE OPERATOR CLASS ecircle_ops
+ DEFAULT FOR TYPE ecircle USING gist AS
+ OPERATOR 13 = ,
+ OPERATOR 21 && (ecircle, epoint),
+ OPERATOR 22 && (ecircle, ebox),
+ OPERATOR 122 &&+ (ecircle, ebox),
+ OPERATOR 23 && (ecircle, ecircle),
+ OPERATOR 24 && (ecircle, ecluster),
+ OPERATOR 124 &&+ (ecircle, ecluster),
+ OPERATOR 31 <-> (ecircle, epoint) FOR ORDER BY float_ops,
+ OPERATOR 32 <-> (ecircle, ebox) FOR ORDER BY float_ops,
+ OPERATOR 33 <-> (ecircle, ecircle) FOR ORDER BY float_ops,
+ OPERATOR 34 <-> (ecircle, ecluster) FOR ORDER BY float_ops,
+ FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
+ FUNCTION 2 pgl_gist_union(internal, internal),
+ FUNCTION 3 pgl_gist_compress_ecircle(internal),
+ FUNCTION 4 pgl_gist_decompress(internal),
+ FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
+ FUNCTION 6 pgl_gist_picksplit(internal, internal),
+ FUNCTION 7 pgl_gist_same(internal, internal, internal),
+ FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
+ STORAGE ekey_area;
+
+CREATE OPERATOR CLASS ecluster_ops
+ DEFAULT FOR TYPE ecluster USING gist AS
+ OPERATOR 21 && (ecluster, epoint),
+ OPERATOR 121 &&+ (ecluster, epoint),
+ OPERATOR 221 @> (ecluster, epoint),
+ OPERATOR 22 && (ecluster, ebox),
+ OPERATOR 122 &&+ (ecluster, ebox),
+ OPERATOR 222 @> (ecluster, ebox),
+ OPERATOR 322 <@ (ecluster, ebox),
+ OPERATOR 23 && (ecluster, ecircle),
+ OPERATOR 123 &&+ (ecluster, ecircle),
+ OPERATOR 24 && (ecluster, ecluster),
+ OPERATOR 124 &&+ (ecluster, ecluster),
+ OPERATOR 224 @> (ecluster, ecluster),
+ OPERATOR 324 <@ (ecluster, ecluster),
+ OPERATOR 31 <-> (ecluster, epoint) FOR ORDER BY float_ops,
+ OPERATOR 32 <-> (ecluster, ebox) FOR ORDER BY float_ops,
+ OPERATOR 33 <-> (ecluster, ecircle) FOR ORDER BY float_ops,
+ OPERATOR 34 <-> (ecluster, ecluster) FOR ORDER BY float_ops,
+ FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
+ FUNCTION 2 pgl_gist_union(internal, internal),
+ FUNCTION 3 pgl_gist_compress_ecluster(internal),
+ FUNCTION 4 pgl_gist_decompress(internal),
+ FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
+ FUNCTION 6 pgl_gist_picksplit(internal, internal),
+ FUNCTION 7 pgl_gist_same(internal, internal, internal),
+ FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
+ STORAGE ekey_area;
+
+
+---------------------
+-- alias functions --
+---------------------
+
+CREATE FUNCTION distance(epoint, epoint)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
+
+CREATE FUNCTION distance(ecluster, epoint)
+ RETURNS float8
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
+
+CREATE FUNCTION distance_within(epoint, epoint, float8)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
+
+CREATE FUNCTION distance_within(ecluster, epoint, float8)
+ RETURNS boolean
+ LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
+
+
+--------------------------------
+-- other data storage formats --
+--------------------------------
+
+CREATE FUNCTION coords_to_epoint(float8, float8, text = 'epoint')
+ RETURNS epoint
+ LANGUAGE plpgsql IMMUTABLE STRICT AS $$
+ DECLARE
+ "result" epoint;
+ BEGIN
+ IF $3 = 'epoint_lonlat' THEN
+ -- avoid dynamic command execution for better performance
+ RETURN epoint($2, $1);
+ END IF;
+ IF $3 = 'epoint' OR $3 = 'epoint_latlon' THEN
+ -- avoid dynamic command execution for better performance
+ RETURN epoint($1, $2);
+ END IF;
+ EXECUTE 'SELECT ' || $3 || '($1, $2)' INTO STRICT "result" USING $1, $2;
+ RETURN "result";
+ END;
+ $$;
+
+CREATE FUNCTION GeoJSON_LinearRing_vertices(jsonb, text = 'epoint_lonlat')
+ RETURNS SETOF jsonb
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT "result" FROM
+ ( SELECT jsonb_array_length($1) - 1 ) AS "lastindex_row" ("lastindex")
+ CROSS JOIN LATERAL jsonb_array_elements(
+ CASE WHEN
+ coords_to_epoint(
+ ($1->0->>0)::float8,
+ ($1->0->>1)::float8,
+ $2
+ ) = coords_to_epoint(
+ ($1->"lastindex"->>0)::float8,
+ ($1->"lastindex"->>1)::float8,
+ $2
+ )
+ THEN
+ $1 - "lastindex"
+ ELSE
+ $1
+ END
+ ) AS "result_row" ("result")
+ $$;
+
+CREATE FUNCTION GeoJSON_to_epoint(jsonb, text = 'epoint_lonlat')
+ RETURNS epoint
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT CASE
+ WHEN $1->>'type' = 'Point' THEN
+ coords_to_epoint(
+ ($1->'coordinates'->>0)::float8,
+ ($1->'coordinates'->>1)::float8,
+ $2
+ )
+ WHEN $1->>'type' = 'Feature' THEN
+ GeoJSON_to_epoint($1->'geometry', $2)
+ ELSE
+ NULL
+ END
+ $$;
+
+CREATE FUNCTION GeoJSON_to_ecluster(jsonb, text = 'epoint_lonlat')
+ RETURNS ecluster
+ LANGUAGE sql IMMUTABLE STRICT AS $$
+ SELECT CASE $1->>'type'
+ WHEN 'Point' THEN
+ coords_to_epoint(
+ ($1->'coordinates'->>0)::float8,
+ ($1->'coordinates'->>1)::float8,
+ $2
+ )::ecluster
+ WHEN 'MultiPoint' THEN
+ ( SELECT ecluster_create_multipoint(array_agg(
+ coords_to_epoint(
+ ("coord"->>0)::float8,
+ ("coord"->>1)::float8,
+ $2
+ )
+ ))
+ FROM jsonb_array_elements($1->'coordinates') AS "coord"
+ )
+ WHEN 'LineString' THEN
+ ( SELECT ecluster_create_path(array_agg(
+ coords_to_epoint(
+ ("coord"->>0)::float8,
+ ("coord"->>1)::float8,
+ $2
+ )
+ ))
+ FROM jsonb_array_elements($1->'coordinates') AS "coord"
+ )
+ WHEN 'MultiLineString' THEN
+ ( SELECT ecluster_concat(array_agg(
+ ( SELECT ecluster_create_path(array_agg(
+ coords_to_epoint(
+ ("coord"->>0)::float8,
+ ("coord"->>1)::float8,
+ $2
+ )
+ ))
+ FROM jsonb_array_elements("coord_array") AS "coord"
+ )
+ ))
+ FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
+ )
+ WHEN 'Polygon' THEN
+ ( SELECT ecluster_concat(array_agg(
+ ( SELECT ecluster_create_polygon(array_agg(
+ coords_to_epoint(
+ ("coord"->>0)::float8,
+ ("coord"->>1)::float8,
+ $2
+ )
+ ))
+ FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
+ )
+ ))
+ FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
+ )
+ WHEN 'MultiPolygon' THEN
+ ( SELECT ecluster_concat(array_agg(
+ ( SELECT ecluster_concat(array_agg(
+ ( SELECT ecluster_create_polygon(array_agg(
+ coords_to_epoint(
+ ("coord"->>0)::float8,
+ ("coord"->>1)::float8,
+ $2
+ )
+ ))
+ FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
+ )
+ ))
+ FROM jsonb_array_elements("coord_array_array") AS "coord_array"
+ )
+ ))
+ FROM jsonb_array_elements($1->'coordinates') AS "coord_array_array"
+ )
+ WHEN 'GeometryCollection' THEN
+ ( SELECT ecluster_concat(array_agg(
+ GeoJSON_to_ecluster("geometry", $2)
+ ))
+ FROM jsonb_array_elements($1->'geometries') AS "geometry"
+ )
+ WHEN 'Feature' THEN
+ GeoJSON_to_ecluster($1->'geometry', $2)
+ WHEN 'FeatureCollection' THEN
+ ( SELECT ecluster_concat(array_agg(
+ GeoJSON_to_ecluster("feature", $2)
+ ))
+ FROM jsonb_array_elements($1->'features') AS "feature"
+ )
+ ELSE
+ NULL
+ END
+ $$;
+
diff -r a90c0d750405 -r 1b9cd45e9e48 latlon-v0008.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/latlon-v0008.c Tue Oct 25 18:44:43 2016 +0200
@@ -0,0 +1,3285 @@
+
+/*-------------*
+ * C prelude *
+ *-------------*/
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "libpq/pqformat.h"
+#include "access/gist.h"
+#include "access/stratnum.h"
+#include "utils/array.h"
+#include
+
+#ifdef PG_MODULE_MAGIC
+PG_MODULE_MAGIC;
+#endif
+
+#if INT_MAX < 2147483647
+#error Expected int type to be at least 32 bit wide
+#endif
+
+
+/*---------------------------------*
+ * distance calculation on earth *
+ * (using WGS-84 spheroid) *
+ *---------------------------------*/
+
+/* WGS-84 spheroid with following parameters:
+ semi-major axis a = 6378137
+ semi-minor axis b = a * (1 - 1/298.257223563)
+ estimated diameter = 2 * (2*a+b)/3
+*/
+#define PGL_SPHEROID_A 6378137.0 /* semi major axis */
+#define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */
+#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
+#define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
+ PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
+ ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
+#define PGL_SUBEPS2 (1.0-PGL_EPS2)
+#define PGL_RADIUS ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0)
+#define PGL_DIAMETER (2.0 * PGL_RADIUS)
+#define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */
+#define PGL_MAXDIST (PGL_RADIUS * M_PI) /* maximum distance */
+#define PGL_FADELIMIT (PGL_MAXDIST / 3.0) /* 1/6 circumference */
+
+/* calculate distance between two points on earth (given in degrees) */
+static inline double pgl_distance(
+ double lat1, double lon1, double lat2, double lon2
+) {
+ float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
+ float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
+ /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
+ /* lon1 = 0 (not used anymore) */
+ lon2 = fabs(lon2-lon1);
+ /* convert to radians (first divide, then multiply) */
+ lat1 = (lat1 / 180.0) * M_PI;
+ lat2 = (lat2 / 180.0) * M_PI;
+ lon2 = (lon2 / 180.0) * M_PI;
+ /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
+ operations (lon2 >= lon1 is already ensured in a previous step) */
+ if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
+ /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
+ of 1.0 */
+ lat1cos = cos(lat1); lat1sin = sin(lat1);
+ lat2cos = cos(lat2); lat2sin = sin(lat2);
+ lon2cos = cos(lon2); lon2sin = sin(lon2);
+ nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
+ nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
+ x1 = nphi1 * lat1cos;
+ z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
+ x2 = nphi2 * lat2cos * lon2cos;
+ y2 = nphi2 * lat2cos * lon2sin;
+ z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
+ /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
+ g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
+ /* convert tunnel distance through scaled ellipsoid to approximated surface
+ distance on original ellipsoid */
+ if (g > 1.0) g = 1.0;
+ s = PGL_DIAMETER * asin(g);
+ /* return result only if small enough to be precise (less than 1/3 of
+ maximum possible distance) */
+ if (s <= PGL_FADELIMIT) return s;
+ /* calculate tunnel distance to antipodal point through scaled ellipsoid */
+ g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
+ /* convert tunnel distance to antipodal point through scaled ellipsoid to
+ approximated surface distance to antipodal point on original ellipsoid */
+ if (g > 1.0) g = 1.0;
+ t = PGL_DIAMETER * asin(g);
+ /* surface distance between original points can now be approximated by
+ substracting antipodal distance from maximum possible distance;
+ return result only if small enough (less than 1/3 of maximum possible
+ distance) */
+ if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
+ /* otherwise crossfade direct and antipodal result to ensure monotonicity */
+ return (
+ (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
+ (s + t - 2*PGL_FADELIMIT)
+ );
+}
+
+/* finite distance that can not be reached on earth */
+#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
+
+
+/*--------------------------------*
+ * simple geographic data types *
+ *--------------------------------*/
+
+/* point on earth given by latitude and longitude in degrees */
+/* (type "epoint" in SQL) */
+typedef struct {
+ double lat; /* between -90 and 90 (both inclusive) */
+ double lon; /* between -180 and 180 (both inclusive) */
+} pgl_point;
+
+/* box delimited by two parallels and two meridians (all in degrees) */
+/* (type "ebox" in SQL) */
+typedef struct {
+ double lat_min; /* between -90 and 90 (both inclusive) */
+ double lat_max; /* between -90 and 90 (both inclusive) */
+ double lon_min; /* between -180 and 180 (both inclusive) */
+ double lon_max; /* between -180 and 180 (both inclusive) */
+ /* if lat_min > lat_max, then box is empty */
+ /* if lon_min > lon_max, then 180th meridian is crossed */
+} pgl_box;
+
+/* circle on earth surface (for radial searches with fixed radius) */
+/* (type "ecircle" in SQL) */
+typedef struct {
+ pgl_point center;
+ double radius; /* positive (including +0 but excluding -0), or -INFINITY */
+ /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
+ zero radius (0) denotes a single point,
+ a finite radius (0 < radius < INFINITY) denotes a filled circle, and
+ a radius of INFINITY is valid and means complete coverage of earth. */
+} pgl_circle;
+
+
+/*----------------------------------*
+ * geographic "cluster" data type *
+ *----------------------------------*/
+
+/* A cluster is a collection of points, paths, outlines, and polygons. If two
+ polygons in a cluster overlap, the area covered by both polygons does not
+ belong to the cluster. This way, a cluster can be used to describe complex
+ shapes like polygons with holes. Outlines are non-filled polygons. Paths are
+ open by default (i.e. the last point in the list is not connected with the
+ first point in the list). Note that each outline or polygon in a cluster
+ must cover a longitude range of less than 180 degrees to avoid ambiguities.
+ Areas which are larger may be split into multiple polygons. */
+
+/* maximum number of points in a cluster */
+/* (limited to avoid integer overflows, e.g. when allocating memory) */
+#define PGL_CLUSTER_MAXPOINTS 16777216
+
+/* types of cluster entries */
+#define PGL_ENTRY_POINT 1 /* a point */
+#define PGL_ENTRY_PATH 2 /* a path from first point to last point */
+#define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */
+#define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */
+
+/* Entries of a cluster are described by two different structs: pgl_newentry
+ and pgl_entry. The first is used only during construction of a cluster, the
+ second is used in all other cases (e.g. when reading clusters from the
+ database, performing operations, etc). */
+
+/* entry for new geographic cluster during construction of that cluster */
+typedef struct {
+ int32_t entrytype;
+ int32_t npoints;
+ pgl_point *points; /* pointer to an array of points (pgl_point) */
+} pgl_newentry;
+
+/* entry of geographic cluster */
+typedef struct {
+ int32_t entrytype; /* type of entry: point, path, outline, polygon */
+ int32_t npoints; /* number of stored points (set to 1 for point entry) */
+ int32_t offset; /* offset of pgl_point array from cluster base address */
+ /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
+} pgl_entry;
+
+/* geographic cluster which is a collection of points, (open) paths, polygons,
+ and outlines (non-filled polygons) */
+typedef struct {
+ char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */
+ int32_t nentries; /* number of stored points */
+ pgl_circle bounding; /* bounding circle */
+ /* Note: bounding circle ensures alignment of pgl_cluster for points */
+ pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */
+} pgl_cluster;
+
+/* macro to determine memory alignment of points */
+/* (needed to store pgl_point array after entries in pgl_cluster) */
+typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
+#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
+
+/* macro to extract a pointer to the array of points of a cluster entry */
+#define PGL_ENTRY_POINTS(cluster, idx) \
+ ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
+
+/* convert pgl_newentry array to pgl_cluster */
+/* NOTE: requires pgl_finalize_cluster to be called to finalize result */
+static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
+ int i; /* index of current entry */
+ int npoints = 0; /* number of points in whole cluster */
+ int entry_npoints; /* number of points in current entry */
+ int points_offset = PGL_POINT_ALIGNMENT * (
+ ( offsetof(pgl_cluster, entries) +
+ nentries * sizeof(pgl_entry) +
+ PGL_POINT_ALIGNMENT - 1
+ ) / PGL_POINT_ALIGNMENT
+ ); /* offset of pgl_point array from base address (considering alignment) */
+ pgl_cluster *cluster; /* new cluster to be returned */
+ /* determine total number of points */
+ for (i=0; ientries[i].entrytype = entries[i].entrytype;
+ cluster->entries[i].npoints = entry_npoints;
+ /* calculate offset (in bytes) of pgl_point array */
+ cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
+ /* copy points */
+ memcpy(
+ PGL_ENTRY_POINTS(cluster, i),
+ entries[i].points,
+ entry_npoints * sizeof(pgl_point)
+ );
+ /* update total number of points processed */
+ npoints += entry_npoints;
+ }
+ /* set number of entries in cluster */
+ cluster->nentries = nentries;
+ /* set PostgreSQL header for variable sized data */
+ SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
+ /* return newly created cluster */
+ return cluster;
+}
+
+
+/*----------------------------------------*
+ * C functions on geographic data types *
+ *----------------------------------------*/
+
+/* round latitude or longitude to 12 digits after decimal point */
+static inline double pgl_round(double val) {
+ return round(val * 1e12) / 1e12;
+}
+
+/* compare two points */
+/* (equality when same point on earth is described, otherwise an arbitrary
+ linear order) */
+static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
+ double lon1, lon2; /* modified longitudes for special cases */
+ /* use latitude as first ordering criterion */
+ if (point1->lat < point2->lat) return -1;
+ if (point1->lat > point2->lat) return 1;
+ /* determine modified longitudes (considering special case of poles and
+ 180th meridian which can be described as W180 or E180) */
+ if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
+ else if (point1->lon == 180) lon1 = -180;
+ else lon1 = point1->lon;
+ if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
+ else if (point2->lon == 180) lon2 = -180;
+ else lon2 = point2->lon;
+ /* use (modified) longitude as secondary ordering criterion */
+ if (lon1 < lon2) return -1;
+ if (lon1 > lon2) return 1;
+ /* no difference found, points are equal */
+ return 0;
+}
+
+/* compare two boxes */
+/* (equality when same box on earth is described, otherwise an arbitrary linear
+ order) */
+static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
+ /* two empty boxes are equal, and an empty box is always considered "less
+ than" a non-empty box */
+ if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
+ if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
+ if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
+ /* use southern border as first ordering criterion */
+ if (box1->lat_min < box2->lat_min) return -1;
+ if (box1->lat_min > box2->lat_min) return 1;
+ /* use northern border as second ordering criterion */
+ if (box1->lat_max < box2->lat_max) return -1;
+ if (box1->lat_max > box2->lat_max) return 1;
+ /* use western border as third ordering criterion */
+ if (box1->lon_min < box2->lon_min) return -1;
+ if (box1->lon_min > box2->lon_min) return 1;
+ /* use eastern border as fourth ordering criterion */
+ if (box1->lon_max < box2->lon_max) return -1;
+ if (box1->lon_max > box2->lon_max) return 1;
+ /* no difference found, boxes are equal */
+ return 0;
+}
+
+/* compare two circles */
+/* (equality when same circle on earth is described, otherwise an arbitrary
+ linear order) */
+static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
+ /* two circles with same infinite radius (positive or negative infinity) are
+ considered equal independently of center point */
+ if (
+ !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
+ circle1->radius == circle2->radius
+ ) return 0;
+ /* use radius as first ordering criterion */
+ if (circle1->radius < circle2->radius) return -1;
+ if (circle1->radius > circle2->radius) return 1;
+ /* use center point as secondary ordering criterion */
+ return pgl_point_cmp(&(circle1->center), &(circle2->center));
+}
+
+/* set box to empty box*/
+static void pgl_box_set_empty(pgl_box *box) {
+ box->lat_min = INFINITY;
+ box->lat_max = -INFINITY;
+ box->lon_min = 0;
+ box->lon_max = 0;
+}
+
+/* check if point is inside a box */
+static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
+ return (
+ point->lat >= box->lat_min && point->lat <= box->lat_max && (
+ (box->lon_min > box->lon_max) ? (
+ /* box crosses 180th meridian */
+ point->lon >= box->lon_min || point->lon <= box->lon_max
+ ) : (
+ /* box does not cross the 180th meridian */
+ point->lon >= box->lon_min && point->lon <= box->lon_max
+ )
+ )
+ );
+}
+
+/* check if two boxes overlap */
+static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
+ return (
+ box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */
+ ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
+ ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
+ (
+ /* check if one and only one box crosses the 180th meridian */
+ ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
+ ((box2->lon_min > box2->lon_max) ? 1 : 0)
+ ) ? (
+ /* exactly one box crosses the 180th meridian */
+ box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
+ box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
+ ) : (
+ /* no box or both boxes cross the 180th meridian */
+ (
+ (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
+ (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
+ ) ||
+ /* handle W180 == E180 */
+ ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
+ ( box2->lon_min == -180 && box1->lon_max == 180 )
+ )
+ )
+ );
+}
+
+/* check unambiguousness of east/west orientation of cluster entries and set
+ bounding circle of cluster */
+static bool pgl_finalize_cluster(pgl_cluster *cluster) {
+ int i, j; /* i: index of entry, j: index of point in entry */
+ int npoints; /* number of points in entry */
+ int total_npoints = 0; /* total number of points in cluster */
+ pgl_point *points; /* points in entry */
+ int lon_dir; /* first point of entry west (-1) or east (+1) */
+ double lon_break = 0; /* antipodal longitude of first point in entry */
+ double lon_min, lon_max; /* covered longitude range of entry */
+ double value; /* temporary variable */
+ /* reset bounding circle center to empty circle at 0/0 coordinates */
+ cluster->bounding.center.lat = 0;
+ cluster->bounding.center.lon = 0;
+ cluster->bounding.radius = -INFINITY;
+ /* if cluster is not empty */
+ if (cluster->nentries != 0) {
+ /* iterate over all cluster entries and ensure they each cover a longitude
+ range less than 180 degrees */
+ for (i=0; inentries; i++) {
+ /* get properties of entry */
+ npoints = cluster->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster, i);
+ /* get longitude of first point of entry */
+ value = points[0].lon;
+ /* initialize lon_min and lon_max with longitude of first point */
+ lon_min = value;
+ lon_max = value;
+ /* determine east/west orientation of first point and calculate antipodal
+ longitude (Note: rounding required here) */
+ if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
+ else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); }
+ else lon_dir = 0;
+ /* iterate over all other points in entry */
+ for (j=1; jlon_break) value = pgl_round(value - 360);
+ else if (lon_dir>0 && value lon_max) lon_max = value;
+ /* return false if 180 degrees or more are covered */
+ if (lon_max - lon_min >= 180) return false;
+ }
+ }
+ /* iterate over all points of all entries and calculate arbitrary center
+ point for bounding circle (best if center point minimizes the radius,
+ but some error is allowed here) */
+ for (i=0; inentries; i++) {
+ /* get properties of entry */
+ npoints = cluster->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster, i);
+ /* check if first entry */
+ if (i==0) {
+ /* get longitude of first point of first entry in whole cluster */
+ value = points[0].lon;
+ /* initialize lon_min and lon_max with longitude of first point of
+ first entry in whole cluster (used to determine if whole cluster
+ covers a longitude range of 180 degrees or more) */
+ lon_min = value;
+ lon_max = value;
+ /* determine east/west orientation of first point and calculate
+ antipodal longitude (Note: rounding not necessary here) */
+ if (value < 0) { lon_dir = -1; lon_break = value + 180; }
+ else if (value > 0) { lon_dir = 1; lon_break = value - 180; }
+ else lon_dir = 0;
+ }
+ /* iterate over all points in entry */
+ for (j=0; j lon_break) value -= 360;
+ else if (lon_dir > 0 && value < lon_break) value += 360;
+ if (value < lon_min) lon_min = value;
+ else if (value > lon_max) lon_max = value;
+ /* set bounding circle to cover whole earth if more than 180 degrees
+ are covered */
+ if (lon_max - lon_min >= 180) {
+ cluster->bounding.center.lat = 0;
+ cluster->bounding.center.lon = 0;
+ cluster->bounding.radius = INFINITY;
+ return true;
+ }
+ /* add point to bounding circle center (for average calculation) */
+ cluster->bounding.center.lat += points[j].lat;
+ cluster->bounding.center.lon += value;
+ }
+ /* count total number of points */
+ total_npoints += npoints;
+ }
+ /* determine average latitude and longitude of cluster */
+ cluster->bounding.center.lat /= total_npoints;
+ cluster->bounding.center.lon /= total_npoints;
+ /* normalize longitude of center of cluster bounding circle */
+ if (cluster->bounding.center.lon < -180) {
+ cluster->bounding.center.lon += 360;
+ }
+ else if (cluster->bounding.center.lon > 180) {
+ cluster->bounding.center.lon -= 360;
+ }
+ /* round bounding circle center (useful if it is used by other functions) */
+ cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
+ cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
+ /* calculate radius of bounding circle */
+ for (i=0; inentries; i++) {
+ npoints = cluster->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster, i);
+ for (j=0; jbounding.center.lat, cluster->bounding.center.lon,
+ points[j].lat, points[j].lon
+ );
+ if (value > cluster->bounding.radius) cluster->bounding.radius = value;
+ }
+ }
+ }
+ /* return true (east/west orientation is unambiguous) */
+ return true;
+}
+
+/* check if point is inside cluster */
+/* (if point is on perimeter, then true is returned if and only if
+ strict == false) */
+static bool pgl_point_in_cluster(
+ pgl_point *point,
+ pgl_cluster *cluster,
+ bool strict
+) {
+ int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
+ int entrytype; /* type of entry */
+ int npoints; /* number of points in entry */
+ pgl_point *points; /* array of points in entry */
+ int lon_dir = 0; /* first vertex west (-1) or east (+1) */
+ double lon_break = 0; /* antipodal longitude of first vertex */
+ double lat0 = point->lat; /* latitude of point */
+ double lon0; /* (adjusted) longitude of point */
+ double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
+ double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
+ double lon; /* longitude of intersection */
+ int counter = 0; /* counter for intersections east of point */
+ /* iterate over all entries */
+ for (i=0; inentries; i++) {
+ /* get type of entry */
+ entrytype = cluster->entries[i].entrytype;
+ /* skip all entries but polygons if perimeters are excluded */
+ if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
+ /* get points of entry */
+ npoints = cluster->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster, i);
+ /* determine east/west orientation of first point of entry and calculate
+ antipodal longitude */
+ lon_break = points[0].lon;
+ if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
+ else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
+ else lon_dir = 0;
+ /* get longitude of point */
+ lon0 = point->lon;
+ /* consider longitude wrap-around for point */
+ if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
+ else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
+ /* iterate over all edges and vertices */
+ for (j=0; j lon_break) lon1 = pgl_round(lon1 - 360);
+ else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
+ }
+ /* get latitude and longitude of next vertex */
+ lat2 = points[k].lat;
+ lon2 = points[k].lon;
+ /* consider longitude wrap-around for next vertex */
+ if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
+ else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
+ /* return if point is on horizontal (west to east) edge of polygon */
+ if (
+ lat0 == lat1 && lat0 == lat2 &&
+ ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
+ ) return !strict;
+ /* check if edge crosses east/west line of point */
+ if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
+ /* calculate longitude of intersection */
+ lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
+ /* return if intersection goes (approximately) through point */
+ if (pgl_round(lon) == lon0) return !strict;
+ /* count intersection if east of point and entry is polygon*/
+ if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
+ }
+ }
+ }
+ /* return true if number of intersections is odd */
+ return counter & 1;
+}
+
+/* check if all points of the second cluster are strictly inside the first
+ cluster */
+static inline bool pgl_all_cluster_points_strictly_in_cluster(
+ pgl_cluster *outer, pgl_cluster *inner
+) {
+ int i, j; /* i: entry, j: point in entry */
+ int npoints; /* number of points in entry */
+ pgl_point *points; /* array of points in entry */
+ /* iterate over all entries of "inner" cluster */
+ for (i=0; inentries; i++) {
+ /* get properties of entry */
+ npoints = inner->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(inner, i);
+ /* iterate over all points in entry of "inner" cluster */
+ for (j=0; jnentries; i++) {
+ /* get properties of entry */
+ npoints = inner->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(inner, i);
+ /* iterate over all points in entry of "inner" cluster */
+ for (j=0; jnentries; i1++) {
+ /* get properties of entry in cluster1 and skip points */
+ npoints1 = cluster1->entries[i1].npoints;
+ if (npoints1 < 2) continue;
+ entrytype1 = cluster1->entries[i1].entrytype;
+ points1 = PGL_ENTRY_POINTS(cluster1, i1);
+ /* determine east/west orientation of first point and calculate antipodal
+ longitude */
+ lon_break1 = points1[0].lon;
+ if (lon_break1 < 0) {
+ lon_dir1 = -1;
+ lon_break1 = pgl_round(lon_break1 + 180);
+ } else if (lon_break1 > 0) {
+ lon_dir1 = 1;
+ lon_break1 = pgl_round(lon_break1 - 180);
+ } else lon_dir1 = 0;
+ /* iterate over all edges and vertices in cluster1 */
+ for (j1=0; j1lon_break1) lon11 = pgl_round(lon11-360);
+ else if (lon_dir1>0 && lon11lon_break1) lon12 = pgl_round(lon12-360);
+ else if (lon_dir1>0 && lon12nentries; i2++) {
+ /* get points and number of points of entry in cluster2 */
+ npoints2 = cluster2->entries[i2].npoints;
+ if (npoints2 < 2) continue;
+ entrytype2 = cluster2->entries[i2].entrytype;
+ points2 = PGL_ENTRY_POINTS(cluster2, i2);
+ /* determine east/west orientation of first point and calculate antipodal
+ longitude */
+ lon_break2 = points2[0].lon;
+ if (lon_break2 < 0) {
+ lon_dir2 = -1;
+ lon_break2 = pgl_round(lon_break2 + 180);
+ } else if (lon_break2 > 0) {
+ lon_dir2 = 1;
+ lon_break2 = pgl_round(lon_break2 - 180);
+ } else lon_dir2 = 0;
+ /* iterate over all edges and vertices in cluster2 */
+ for (j2=0; j2lon_break2) lon21 = pgl_round(lon21-360);
+ else if (lon_dir2>0 && lon21lon_break2) lon22 = pgl_round(lon22-360);
+ else if (lon_dir2>0 && lon22 360) {
+ lon21 = pgl_round(lon21 - 360);
+ lon22 = pgl_round(lon22 - 360);
+ } else if (wrapvalue < -360) {
+ lon21 = pgl_round(lon21 + 360);
+ lon22 = pgl_round(lon22 + 360);
+ }
+ /* return true if segments overlap */
+ if (
+ pgl_lseg_crosses_line(
+ lat11, lon11, lat12, lon12,
+ lat21, lon21, lat22, lon22
+ ) && pgl_lseg_crosses_line(
+ lat21, lon21, lat22, lon22,
+ lat11, lon11, lat12, lon12
+ )
+ ) {
+ return true;
+ }
+ }
+ }
+ }
+ }
+ /* otherwise return false */
+ return false;
+}
+
+/* check if second cluster is completely contained in first cluster */
+static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
+ if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
+ if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
+ if (pgl_outlines_overlap(outer, inner)) return false;
+ return true;
+}
+
+/* check if two clusters overlap */
+static bool pgl_clusters_overlap(
+ pgl_cluster *cluster1, pgl_cluster *cluster2
+) {
+ if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
+ if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
+ if (pgl_outlines_overlap(cluster1, cluster2)) return true;
+ return false;
+}
+
+
+/* calculate (approximate) distance between point and cluster */
+static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
+ double comp; /* square of compression of meridians */
+ int i, j, k; /* i: entry, j: point in entry, k: next point in entry */
+ int entrytype; /* type of entry */
+ int npoints; /* number of points in entry */
+ pgl_point *points; /* array of points in entry */
+ int lon_dir = 0; /* first vertex west (-1) or east (+1) */
+ double lon_break = 0; /* antipodal longitude of first vertex */
+ double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */
+ double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */
+ double lat0 = point->lat; /* latitude of point */
+ double lon0; /* (adjusted) longitude of point */
+ double lat1, lon1; /* latitude and (adjusted) longitude of vertex */
+ double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */
+ double s; /* scalar for vector calculations */
+ double dist; /* distance calculated in one step */
+ double min_dist = INFINITY; /* minimum distance */
+ /* distance is zero if point is contained in cluster */
+ if (pgl_point_in_cluster(point, cluster, false)) return 0;
+ /* calculate approximate square compression of meridians */
+ comp = cos((lat0 / 180.0) * M_PI);
+ comp *= comp;
+ /* calculate exact square compression of meridians */
+ comp *= (
+ (1.0 - PGL_EPS2 * (1.0-comp)) *
+ (1.0 - PGL_EPS2 * (1.0-comp)) /
+ (PGL_SUBEPS2 * PGL_SUBEPS2)
+ );
+ /* iterate over all entries */
+ for (i=0; inentries; i++) {
+ /* get properties of entry */
+ entrytype = cluster->entries[i].entrytype;
+ npoints = cluster->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster, i);
+ /* determine east/west orientation of first point of entry and calculate
+ antipodal longitude */
+ lon_break = points[0].lon;
+ if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
+ else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
+ else lon_dir = 0;
+ /* determine covered longitude range */
+ for (j=0; j lon_break) lon1 -= 360;
+ else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
+ /* update minimum and maximum longitude of polygon */
+ if (j == 0 || lon1 < lon_min) lon_min = lon1;
+ if (j == 0 || lon1 > lon_max) lon_max = lon1;
+ }
+ /* adjust longitude wrap-around according to full longitude range */
+ lon_break = (lon_max + lon_min) / 2;
+ if (lon_break < 0) { lon_dir = -1; lon_break += 180; }
+ else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; }
+ /* get longitude of point */
+ lon0 = point->lon;
+ /* consider longitude wrap-around for point */
+ if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
+ else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
+ /* iterate over all edges and vertices */
+ for (j=0; j lon_break) lon1 -= 360;
+ else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
+ }
+ /* calculate distance to vertex */
+ dist = pgl_distance(lat0, lon0, lat1, lon1);
+ /* store calculated distance if smallest */
+ if (dist < min_dist) min_dist = dist;
+ /* calculate index of next vertex */
+ k = (j+1) % npoints;
+ /* skip last edge unless entry is (closed) outline or polygon */
+ if (
+ k == 0 &&
+ entrytype != PGL_ENTRY_OUTLINE &&
+ entrytype != PGL_ENTRY_POLYGON
+ ) continue;
+ /* get latitude and longitude of next vertex */
+ lat2 = points[k].lat;
+ lon2 = points[k].lon;
+ /* consider longitude wrap-around for next vertex */
+ if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
+ else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
+ /* go to next vertex and edge if edge is degenerated */
+ if (lat1 == lat2 && lon1 == lon2) continue;
+ /* otherwise test if point can be projected onto edge of polygon */
+ s = (
+ ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
+ ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
+ );
+ /* go to next vertex and edge if point cannot be projected */
+ if (!(s > 0 && s < 1)) continue;
+ /* calculate distance from original point to projected point */
+ dist = pgl_distance(
+ lat0, lon0,
+ lat1 + s * (lat2-lat1),
+ lon1 + s * (lon2-lon1)
+ );
+ /* store calculated distance if smallest */
+ if (dist < min_dist) min_dist = dist;
+ }
+ }
+ /* return minimum distance */
+ return min_dist;
+}
+
+/* calculate (approximate) distance between two clusters */
+static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
+ int i, j; /* i: entry, j: point in entry */
+ int npoints; /* number of points in entry */
+ pgl_point *points; /* array of points in entry */
+ double dist; /* distance calculated in one step */
+ double min_dist = INFINITY; /* minimum distance */
+ /* consider distance from each point in one cluster to the whole other */
+ for (i=0; inentries; i++) {
+ npoints = cluster1->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster1, i);
+ for (j=0; jnentries; i++) {
+ npoints = cluster2->entries[i].npoints;
+ points = PGL_ENTRY_POINTS(cluster2, i);
+ for (j=0; jlat_min > box->lat_max) return INFINITY;
+ /* return zero if point is inside box */
+ if (pgl_point_in_box(point, box)) return 0;
+ /* calculate delta longitude */
+ dlon = box->lon_max - box->lon_min;
+ if (dlon < 0) dlon += 360; /* 180th meridian crossed */
+ /* if delta longitude is greater than 150 degrees, perform safe fall-back */
+ if (dlon > 150) return 0;
+ /* calculate lower limit for distance (formula below requires dlon <= 150) */
+ /* TODO: provide better estimation function to improve performance */
+ distance = (
+ (1.0-1e-14) * pgl_distance(
+ point->lat,
+ point->lon,
+ (box->lat_min + box->lat_max) / 2,
+ box->lon_min + dlon/2
+ ) - pgl_distance(
+ box->lat_min, box->lon_min,
+ box->lat_max, box->lon_max
+ )
+ );
+ /* truncate negative results to zero */
+ if (distance <= 0) distance = 0;
+ /* return result */
+ return distance;
+}
+
+
+/*-------------------------------*
+ * Monte Carlo based functions *
+ *-------------------------------*/
+
+/* half of (spherical) earth's surface area */
+#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI)
+
+/* golden angle in radians */
+#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0))
+
+/* create a list of sample points covering a bounding circle
+ and return covered area */
+static double pgl_sample_points(
+ pgl_point *center, /* center of bounding circle */
+ double radius, /* radius of bounding circle */
+ int samples, /* number of sample points */
+ pgl_point *result /* pointer to result array */
+) {
+ double double_share = 2.0; /* double of covered share of earth's surface */
+ double double_share_div_samples; /* double_share divided by sample count */
+ int i;
+ double t; /* parameter of spiral laid on (spherical) earth's surface */
+ double x, y, z; /* normalized coordinates of point on non-rotated spiral */
+ double sin_phi; /* sine of sph. coordinate of point of non-rotated spiral */
+ double lambda; /* other sph. coordinate of point of non-rotated spiral */
+ double rot = (0.5 - center->lat / 180.0) * M_PI; /* needed rot. (in rad) */
+ double cos_rot = cos(rot); /* cosine of rotation by latitude */
+ double sin_rot = sin(rot); /* sine of rotation by latitude */
+ double x_rot, z_rot; /* normalized coordinates of point on rotated spiral */
+ double center_lon = center->lon; /* second rotation in degree */
+ /* add safety margin to bounding circle because of spherical approximation */
+ radius *= PGL_SPHEROID_A / PGL_RADIUS;
+ /* if whole earth is covered, use initialized value, otherwise calculate
+ share of covered area (multiplied by 2) */
+ if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS);
+ /* divide double_share by sample count for later calculations */
+ double_share_div_samples = double_share / samples;
+ /* generate sample points */
+ for (i=0; i PGL_CLUSTER_MAXPOINTS) {
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg(
+ "too many sample points for Monte Carlo method (maximum %i)",
+ PGL_CLUSTER_MAXPOINTS
+ )
+ ));
+ }
+ if (!(cluster->bounding.radius > 0)) return 0;
+ for (i=0; ; i++) {
+ if (i == cluster->nentries) return 0;
+ if (cluster->entries[i].entrytype == PGL_ENTRY_POLYGON) break;
+ }
+ points = palloc(samples * sizeof(pgl_point));
+ area = pgl_sample_points(
+ &cluster->bounding.center,
+ cluster->bounding.radius,
+ samples,
+ points
+ );
+ for (i=0; i PGL_CLUSTER_MAXPOINTS) {
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg(
+ "too many sample points for Monte Carlo method (maximum %i)",
+ PGL_CLUSTER_MAXPOINTS
+ )
+ ));
+ }
+ /* calculate shortest distance from point to cluster */
+ distance = pgl_point_cluster_distance(point, cluster);
+ /* if cluster consists of a single point or has no bounding circle with
+ positive radius, simply return distance */
+ if (
+ (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) ||
+ !(cluster->bounding.radius > 0)
+ ) return distance;
+ /* if cluster consists of two points which are twice as far apart, return
+ distance between point and cluster multiplied by square root of two */
+ if (
+ cluster->nentries == 2 &&
+ cluster->entries[0].entrytype == PGL_ENTRY_POINT &&
+ cluster->entries[1].entrytype == PGL_ENTRY_POINT &&
+ pgl_distance(
+ PGL_ENTRY_POINTS(cluster, 0)[0].lat,
+ PGL_ENTRY_POINTS(cluster, 0)[0].lon,
+ PGL_ENTRY_POINTS(cluster, 1)[0].lat,
+ PGL_ENTRY_POINTS(cluster, 1)[0].lon
+ ) >= 2.0 * distance
+ ) {
+ return distance * M_SQRT2;
+ }
+ /* otherwise create sample points for Monte Carlo method and determine area
+ covered by sample points */
+ points = palloc(samples * sizeof(pgl_point));
+ area = pgl_sample_points(
+ &cluster->bounding.center,
+ cluster->bounding.radius + distance, /* pad bounding circle by distance */
+ samples,
+ points
+ );
+ /* perform Monte Carlo method */
+ if (distance > 0) {
+ /* point is outside cluster */
+ for (i=0; i distance) return result;
+ /* otherwise return distance between point and cluster */
+ else return distance;
+}
+
+
+/*-------------------------------------------------*
+ * geographic index based on space-filling curve *
+ *-------------------------------------------------*/
+
+/* number of bytes used for geographic (center) position in keys */
+#define PGL_KEY_LATLON_BYTELEN 7
+
+/* maximum reference value for logarithmic size of geographic objects */
+#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */
+
+/* pointer to index key (either pgl_pointkey or pgl_areakey) */
+typedef unsigned char *pgl_keyptr;
+
+/* index key for points (objects with zero area) on the spheroid */
+/* bit 0..55: interspersed bits of latitude and longitude,
+ bit 56..57: always zero,
+ bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
+typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
+
+/* index key for geographic objects on spheroid with area greater than zero */
+/* bit 0..55: interspersed bits of latitude and longitude of center point,
+ bit 56: always set to 1,
+ bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
+ bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
+ PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
+ = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
+ (with interspersed bits = 0 and node depth = 0) for keys which
+ cover both empty and non-empty objects */
+
+typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
+
+/* helper macros for reading/writing index keys */
+#define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN
+#define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1)
+#define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8)
+#define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1)
+#define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1)
+#define PGL_AREAKEY_TYPEMASK 0x80
+#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
+#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
+ ( PGL_KEY_LATLONBIT(key1, n) ^ \
+ PGL_KEY_LATLONBIT(key2, n) )
+#define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
+ PGL_AREAKEY_TYPEMASK)
+#define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
+ (PGL_AREAKEY_TYPEMASK-1))
+#define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET])
+#define PGL_KEY_OBJSIZE_EMPTY 126
+#define PGL_KEY_OBJSIZE_UNIVERSAL 127
+#define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \
+ (key)[PGL_KEY_OBJSIZE_OFFSET] == \
+ PGL_KEY_OBJSIZE_EMPTY )
+#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
+ (key)[PGL_KEY_OBJSIZE_OFFSET] == \
+ PGL_KEY_OBJSIZE_UNIVERSAL )
+
+/* set area key to match empty objects only */
+static void pgl_key_set_empty(pgl_keyptr key) {
+ memset(key, 0, sizeof(pgl_areakey));
+ /* Note: setting node depth to maximum is required for picksplit function */
+ key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
+ key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
+}
+
+/* set area key to match any object (including empty objects) */
+static void pgl_key_set_universal(pgl_keyptr key) {
+ memset(key, 0, sizeof(pgl_areakey));
+ key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
+ key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
+}
+
+/* convert a point on earth into a max-depth key to be used in index */
+static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
+ double lat = point->lat;
+ double lon = point->lon;
+ int i;
+ /* clear latitude and longitude bits */
+ memset(key, 0, PGL_KEY_LATLON_BYTELEN);
+ /* set node depth to maximum and type bit to zero */
+ key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
+ /* iterate over all latitude/longitude bit pairs */
+ for (i=0; i= 0) {
+ key[i/4] |= 0x80 >> (2*(i%4));
+ lat *= 2; lat -= 90;
+ } else {
+ lat *= 2; lat += 90;
+ }
+ /* determine longitude bit */
+ if (lon >= 0) {
+ key[i/4] |= 0x80 >> (2*(i%4)+1);
+ lon *= 2; lon -= 180;
+ } else {
+ lon *= 2; lon += 180;
+ }
+ }
+}
+
+/* convert a circle on earth into a max-depth key to be used in an index */
+static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
+ /* handle special case of empty circle */
+ if (circle->radius < 0) {
+ pgl_key_set_empty(key);
+ return;
+ }
+ /* perform same action as for point keys */
+ pgl_point_to_key(&(circle->center), key);
+ /* but overwrite type and node depth to fit area index key */
+ key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
+ /* check if radius is greater than (or equal to) reference size */
+ /* (treat equal values as greater values for numerical safety) */
+ if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
+ /* if yes, set logarithmic size to zero */
+ key[PGL_KEY_OBJSIZE_OFFSET] = 0;
+ } else {
+ /* otherwise, determine logarithmic size iteratively */
+ /* (one step is equivalent to a factor of sqrt(2)) */
+ double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
+ int objsize = 1;
+ while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
+ /* stop when radius is greater than (or equal to) adjusted reference */
+ /* (treat equal values as greater values for numerical safety) */
+ if (circle->radius >= reference) break;
+ reference /= M_SQRT2;
+ objsize++;
+ }
+ /* set logarithmic size to determined value */
+ key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
+ }
+}
+
+/* check if one key is subkey of another key or vice versa */
+static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
+ int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */
+ /* determine smallest depth */
+ int depth1 = PGL_KEY_NODEDEPTH(key1);
+ int depth2 = PGL_KEY_NODEDEPTH(key2);
+ int depth = (depth1 < depth2) ? depth1 : depth2;
+ /* check if keys are area keys (assuming that both keys have same type) */
+ if (PGL_KEY_IS_AREAKEY(key1)) {
+ int j = 0; /* bit offset for logarithmic object size bits */
+ int k = 0; /* bit offset for latitude and longitude */
+ /* fetch logarithmic object size information */
+ int objsize1 = PGL_KEY_OBJSIZE(key1);
+ int objsize2 = PGL_KEY_OBJSIZE(key2);
+ /* handle special cases for empty objects (universal and empty keys) */
+ if (
+ objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
+ objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
+ ) return true;
+ if (
+ objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
+ objsize2 == PGL_KEY_OBJSIZE_EMPTY
+ ) return objsize1 == objsize2;
+ /* iterate through key bits */
+ for (i=0; i j) ||
+ (objsize2 <= j && objsize1 > j)
+ ) {
+ /* bit differs, therefore keys are in separate branches */
+ return false;
+ }
+ /* increase bit counter for object size bits */
+ j++;
+ }
+ /* all other bits describe latitude and longitude */
+ else {
+ /* check if bit differs in both keys */
+ if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
+ /* bit differs, therefore keys are in separate branches */
+ return false;
+ }
+ /* increase bit counter for latitude/longitude bits */
+ k++;
+ }
+ }
+ }
+ /* if not, keys are point keys */
+ else {
+ /* iterate through key bits */
+ for (i=0; i PGL_AREAKEY_MAXOBJSIZE ||
+ objsize2 > PGL_AREAKEY_MAXOBJSIZE
+ ) {
+ if (
+ objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
+ objsize2 == PGL_KEY_OBJSIZE_EMPTY
+ ) pgl_key_set_empty(dst);
+ else pgl_key_set_universal(dst);
+ return;
+ }
+ /* iterate through key bits */
+ for (i=0; i= j && objsize2 >= j) {
+ /* set objsize in destination buffer to indicate that size bit is
+ unset in destination buffer at the current bit position */
+ dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
+ }
+ /* break if object size bit is set in one key only */
+ else if (objsize1 >= j || objsize2 >= j) break;
+ }
+ /* all other bits describe latitude and longitude */
+ else {
+ /* break if bit differs in both keys */
+ if (PGL_KEY_LATLONBIT(dst, k)) {
+ if (!PGL_KEY_LATLONBIT(src, k)) break;
+ /* but set bit in destination buffer if bit is set in both keys */
+ dstbuf[k/8] |= 0x80 >> (k%8);
+ } else if (PGL_KEY_LATLONBIT(src, k)) break;
+ /* increase bit counter for latitude/longitude bits */
+ k++;
+ }
+ }
+ /* set common node depth and type bit (type bit = 1) */
+ dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
+ /* copy contents of destination buffer to first key */
+ memcpy(dst, dstbuf, sizeof(pgl_areakey));
+ }
+ /* if not, keys are point keys */
+ else {
+ pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */
+ /* iterate through key bits */
+ for (i=0; i> (i%8);
+ } else if (PGL_KEY_LATLONBIT(src, i)) break;
+ }
+ /* set common node depth (type bit = 0) */
+ dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
+ /* copy contents of destination buffer to first key */
+ memcpy(dst, dstbuf, sizeof(pgl_pointkey));
+ }
+}
+
+/* determine center(!) boundaries and radius estimation of index key */
+static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
+ int i;
+ /* determine node depth */
+ int depth = PGL_KEY_NODEDEPTH(key);
+ /* center point of possible result */
+ double lat = 0;
+ double lon = 0;
+ /* maximum distance of real center point from key center */
+ double dlat = 90;
+ double dlon = 180;
+ /* maximum radius of contained objects */
+ double radius = 0; /* always return zero for point index keys */
+ /* check if key is area key */
+ if (PGL_KEY_IS_AREAKEY(key)) {
+ /* get logarithmic object size */
+ int objsize = PGL_KEY_OBJSIZE(key);
+ /* handle special cases for empty objects (universal and empty keys) */
+ if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
+ pgl_box_set_empty(box);
+ return 0;
+ } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
+ box->lat_min = -90;
+ box->lat_max = 90;
+ box->lon_min = -180;
+ box->lon_max = 180;
+ return 0; /* any value >= 0 would do */
+ }
+ /* calculate maximum possible radius of objects covered by the given key */
+ if (objsize == 0) radius = INFINITY;
+ else {
+ radius = PGL_AREAKEY_REFOBJSIZE;
+ while (--objsize) radius /= M_SQRT2;
+ }
+ /* iterate over latitude and longitude bits in key */
+ /* (every second bit is a latitude or longitude bit) */
+ for (i=0; ilat_min = lat - dlat;
+ box->lat_max = lat + dlat;
+ box->lon_min = lon - dlon;
+ box->lon_max = lon + dlon;
+ /* return radius (as a function return value) */
+ return radius;
+}
+
+/* estimator function for distance between point and index key */
+/* always returns a smaller value than actually correct or zero */
+static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
+ pgl_box box; /* center(!) bounding box of area index key */
+ /* calculate center(!) bounding box and maximum radius of objects covered
+ by area index key (radius is zero for point index keys) */
+ double distance = pgl_key_to_box(key, &box);
+ /* calculate estimated distance between bounding box of center point of
+ indexed object and point passed as second argument, then substract maximum
+ radius of objects covered by index key */
+ distance = pgl_estimate_point_box_distance(point, &box) - distance;
+ /* truncate negative results to zero */
+ if (distance <= 0) distance = 0;
+ /* return result */
+ return distance;
+}
+
+
+/*---------------------------------*
+ * helper functions for text I/O *
+ *---------------------------------*/
+
+#define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */
+
+/* convert floating point number to string (round-trip safe) */
+static void pgl_print_float(char *buf, double flt) {
+ /* check if number is integral */
+ if (trunc(flt) == flt) {
+ /* for integral floats use maximum precision */
+ snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
+ } else {
+ /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
+ snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
+ if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
+ if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
+ }
+}
+
+/* convert latitude floating point number (in degrees) to string */
+static void pgl_print_lat(char *buf, double lat) {
+ if (signbit(lat)) {
+ /* treat negative latitudes (including -0) as south */
+ snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
+ } else {
+ /* treat positive latitudes (including +0) as north */
+ snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
+ }
+}
+
+/* convert longitude floating point number (in degrees) to string */
+static void pgl_print_lon(char *buf, double lon) {
+ if (signbit(lon)) {
+ /* treat negative longitudes (including -0) as west */
+ snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
+ } else {
+ /* treat positive longitudes (including +0) as east */
+ snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
+ }
+}
+
+/* bit masks used as return value of pgl_scan() function */
+#define PGL_SCAN_NONE 0 /* no value has been parsed */
+#define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */
+#define PGL_SCAN_LON (1<<1) /* longitude has been parsed */
+#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */
+
+/* parse a coordinate (can be latitude or longitude) */
+static int pgl_scan(char **str, double *lat, double *lon) {
+ double val;
+ int len;
+ if (
+ sscanf(*str, " N %lf %n", &val, &len) ||
+ sscanf(*str, " n %lf %n", &val, &len)
+ ) {
+ *str += len; *lat = val; return PGL_SCAN_LAT;
+ }
+ if (
+ sscanf(*str, " S %lf %n", &val, &len) ||
+ sscanf(*str, " s %lf %n", &val, &len)
+ ) {
+ *str += len; *lat = -val; return PGL_SCAN_LAT;
+ }
+ if (
+ sscanf(*str, " E %lf %n", &val, &len) ||
+ sscanf(*str, " e %lf %n", &val, &len)
+ ) {
+ *str += len; *lon = val; return PGL_SCAN_LON;
+ }
+ if (
+ sscanf(*str, " W %lf %n", &val, &len) ||
+ sscanf(*str, " w %lf %n", &val, &len)
+ ) {
+ *str += len; *lon = -val; return PGL_SCAN_LON;
+ }
+ return PGL_SCAN_NONE;
+}
+
+
+/*-----------------*
+ * SQL functions *
+ *-----------------*/
+
+/* Note: These function names use "epoint", "ebox", etc. notation here instead
+ of "point", "box", etc. in order to distinguish them from any previously
+ defined functions. */
+
+/* function needed for dummy types and/or not implemented features */
+PG_FUNCTION_INFO_V1(pgl_notimpl);
+Datum pgl_notimpl(PG_FUNCTION_ARGS) {
+ ereport(ERROR, (errmsg("not implemented by pgLatLon")));
+}
+
+/* set point to latitude and longitude (including checks) */
+static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
+ /* reject infinite or NaN values */
+ if (!isfinite(lat) || !isfinite(lon)) {
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg("epoint requires finite coordinates")
+ ));
+ }
+ /* check latitude bounds */
+ if (lat < -90) {
+ ereport(WARNING, (errmsg("latitude exceeds south pole")));
+ lat = -90;
+ } else if (lat > 90) {
+ ereport(WARNING, (errmsg("latitude exceeds north pole")));
+ lat = 90;
+ }
+ /* check longitude bounds */
+ if (lon < -180) {
+ ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
+ lon += 360 - trunc(lon / 360) * 360;
+ } else if (lon > 180) {
+ ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
+ lon -= 360 + trunc(lon / 360) * 360;
+ }
+ /* store rounded latitude/longitude values for round-trip safety */
+ point->lat = pgl_round(lat);
+ point->lon = pgl_round(lon);
+}
+
+/* create point ("epoint" in SQL) from latitude and longitude */
+PG_FUNCTION_INFO_V1(pgl_create_epoint);
+Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
+ pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
+ PG_RETURN_POINTER(point);
+}
+
+/* parse point ("epoint" in SQL) */
+/* format: '[NS] [EW]' */
+PG_FUNCTION_INFO_V1(pgl_epoint_in);
+Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
+ char *str = PG_GETARG_CSTRING(0); /* input string */
+ char *strptr = str; /* current position within string */
+ int done = 0; /* bit mask storing if latitude or longitude was read */
+ double lat, lon; /* parsed values as double precision floats */
+ pgl_point *point; /* return value (to be palloc'ed) */
+ /* parse two floats (each latitude or longitude) separated by white-space */
+ done |= pgl_scan(&strptr, &lat, &lon);
+ if (strptr != str && isspace(strptr[-1])) {
+ done |= pgl_scan(&strptr, &lat, &lon);
+ }
+ /* require end of string, and latitude and longitude parsed successfully */
+ if (strptr[0] || done != PGL_SCAN_LATLON) {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+ errmsg("invalid input syntax for type epoint: \"%s\"", str)
+ ));
+ }
+ /* allocate memory for result */
+ point = (pgl_point *)palloc(sizeof(pgl_point));
+ /* set latitude and longitude (and perform checks) */
+ pgl_epoint_set_latlon(point, lat, lon);
+ /* return result */
+ PG_RETURN_POINTER(point);
+}
+
+/* create box ("ebox" in SQL) that is empty */
+PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
+Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
+ pgl_box_set_empty(box);
+ PG_RETURN_POINTER(box);
+}
+
+/* set box to given boundaries (including checks) */
+static void pgl_ebox_set_boundaries(
+ pgl_box *box,
+ double lat_min, double lat_max, double lon_min, double lon_max
+) {
+ /* if minimum latitude is greater than maximum latitude, return empty box */
+ if (lat_min > lat_max) {
+ pgl_box_set_empty(box);
+ return;
+ }
+ /* otherwise reject infinite or NaN values */
+ if (
+ !isfinite(lat_min) || !isfinite(lat_max) ||
+ !isfinite(lon_min) || !isfinite(lon_max)
+ ) {
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg("ebox requires finite coordinates")
+ ));
+ }
+ /* check latitude bounds */
+ if (lat_max < -90) {
+ ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
+ lat_max = -90;
+ } else if (lat_max > 90) {
+ ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
+ lat_max = 90;
+ }
+ if (lat_min < -90) {
+ ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
+ lat_min = -90;
+ } else if (lat_min > 90) {
+ ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
+ lat_min = 90;
+ }
+ /* check if all longitudes are included */
+ if (lon_max - lon_min >= 360) {
+ if (lon_max - lon_min > 360) ereport(WARNING, (
+ errmsg("longitude coverage greater than 360 degrees")
+ ));
+ lon_min = -180;
+ lon_max = 180;
+ } else {
+ /* normalize longitude bounds */
+ if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
+ else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360;
+ if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
+ else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360;
+ }
+ /* store rounded latitude/longitude values for round-trip safety */
+ box->lat_min = pgl_round(lat_min);
+ box->lat_max = pgl_round(lat_max);
+ box->lon_min = pgl_round(lon_min);
+ box->lon_max = pgl_round(lon_max);
+ /* ensure that rounding does not change orientation */
+ if (lon_min > lon_max && box->lon_min == box->lon_max) {
+ box->lon_min = -180;
+ box->lon_max = 180;
+ }
+}
+
+/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
+PG_FUNCTION_INFO_V1(pgl_create_ebox);
+Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
+ pgl_ebox_set_boundaries(
+ box,
+ PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
+ PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
+ );
+ PG_RETURN_POINTER(box);
+}
+
+/* create box ("ebox" in SQL) from two points ("epoint"s) */
+/* (can not be used to cover a longitude range of more than 120 degrees) */
+PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
+Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
+ pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
+ pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
+ double lat_min, lat_max, lon_min, lon_max;
+ double dlon; /* longitude range (delta longitude) */
+ /* order latitude and longitude boundaries */
+ if (point2->lat < point1->lat) {
+ lat_min = point2->lat;
+ lat_max = point1->lat;
+ } else {
+ lat_min = point1->lat;
+ lat_max = point2->lat;
+ }
+ if (point2->lon < point1->lon) {
+ lon_min = point2->lon;
+ lon_max = point1->lon;
+ } else {
+ lon_min = point1->lon;
+ lon_max = point2->lon;
+ }
+ /* calculate longitude range (round to avoid floating point errors) */
+ dlon = pgl_round(lon_max - lon_min);
+ /* determine east-west direction */
+ if (dlon >= 240) {
+ /* assume that 180th meridian is crossed and swap min/max longitude */
+ double swap = lon_min; lon_min = lon_max; lon_max = swap;
+ } else if (dlon > 120) {
+ /* unclear orientation since delta longitude > 120 */
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg("can not determine east/west orientation for ebox")
+ ));
+ }
+ /* use boundaries to setup box (and perform checks) */
+ pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
+ /* return result */
+ PG_RETURN_POINTER(box);
+}
+
+/* parse box ("ebox" in SQL) */
+/* format: '[NS] [EW] [NS] [EW]'
+ or: '[NS] [NS] [EW] [EW]' */
+PG_FUNCTION_INFO_V1(pgl_ebox_in);
+Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
+ char *str = PG_GETARG_CSTRING(0); /* input string */
+ char *str_lower; /* lower case version of input string */
+ char *strptr; /* current position within string */
+ int valid; /* number of valid chars */
+ int done; /* specifies if latitude or longitude was read */
+ double val; /* temporary variable */
+ int lat_count = 0; /* count of latitude values parsed */
+ int lon_count = 0; /* count of longitufde values parsed */
+ double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */
+ pgl_box *box; /* return value (to be palloc'ed) */
+ /* lowercase input */
+ str_lower = psprintf("%s", str);
+ for (strptr=str_lower; *strptr; strptr++) {
+ if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
+ }
+ /* reset reading position to start of (lowercase) string */
+ strptr = str_lower;
+ /* check if empty box */
+ valid = 0;
+ sscanf(strptr, " empty %n", &valid);
+ if (valid && strptr[valid] == 0) {
+ /* allocate and return empty box */
+ box = (pgl_box *)palloc(sizeof(pgl_box));
+ pgl_box_set_empty(box);
+ PG_RETURN_POINTER(box);
+ }
+ /* demand four blocks separated by whitespace */
+ valid = 0;
+ sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
+ /* if four blocks separated by whitespace exist, parse those blocks */
+ if (strptr[valid] == 0) while (strptr[0]) {
+ /* parse either latitude or longitude (whichever found in input string) */
+ done = pgl_scan(&strptr, &val, &val);
+ /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
+ if (done == PGL_SCAN_LAT) {
+ if (!lat_count) lat_min = val; else lat_max = val;
+ lat_count++;
+ } else if (done == PGL_SCAN_LON) {
+ if (!lon_count) lon_min = val; else lon_max = val;
+ lon_count++;
+ } else {
+ break;
+ }
+ }
+ /* require end of string, and two latitude and two longitude values */
+ if (strptr[0] || lat_count != 2 || lon_count != 2) {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+ errmsg("invalid input syntax for type ebox: \"%s\"", str)
+ ));
+ }
+ /* free lower case string */
+ pfree(str_lower);
+ /* order boundaries (maximum greater than minimum) */
+ if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
+ if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
+ /* allocate memory for result */
+ box = (pgl_box *)palloc(sizeof(pgl_box));
+ /* set boundaries (and perform checks) */
+ pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
+ /* return result */
+ PG_RETURN_POINTER(box);
+}
+
+/* set circle to given latitude, longitude, and radius (including checks) */
+static void pgl_ecircle_set_latlon_radius(
+ pgl_circle *circle, double lat, double lon, double radius
+) {
+ /* set center point (including checks) */
+ pgl_epoint_set_latlon(&(circle->center), lat, lon);
+ /* handle non-positive radius */
+ if (isnan(radius)) {
+ ereport(ERROR, (
+ errcode(ERRCODE_DATA_EXCEPTION),
+ errmsg("invalid radius for ecircle")
+ ));
+ }
+ if (radius == 0) radius = 0; /* avoids -0 */
+ else if (radius < 0) {
+ if (isfinite(radius)) {
+ ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
+ }
+ radius = -INFINITY;
+ }
+ /* store radius (round-trip safety is ensured by pgl_print_float) */
+ circle->radius = radius;
+}
+
+/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
+PG_FUNCTION_INFO_V1(pgl_create_ecircle);
+Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
+ pgl_ecircle_set_latlon_radius(
+ circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
+ );
+ PG_RETURN_POINTER(circle);
+}
+
+/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
+PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
+Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ double radius = PG_GETARG_FLOAT8(1);
+ pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
+ /* set latitude, longitude, radius (and perform checks) */
+ pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
+ /* return result */
+ PG_RETURN_POINTER(circle);
+}
+
+/* parse circle ("ecircle" in SQL) */
+/* format: '[NS] [EW] ' */
+PG_FUNCTION_INFO_V1(pgl_ecircle_in);
+Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
+ char *str = PG_GETARG_CSTRING(0); /* input string */
+ char *strptr = str; /* current position within string */
+ double lat, lon, radius; /* parsed values as double precision flaots */
+ int valid = 0; /* number of valid chars */
+ int done = 0; /* stores if latitude and/or longitude was read */
+ pgl_circle *circle; /* return value (to be palloc'ed) */
+ /* demand three blocks separated by whitespace */
+ sscanf(strptr, " %*s %*s %*s %n", &valid);
+ /* if three blocks separated by whitespace exist, parse those blocks */
+ if (strptr[valid] == 0) {
+ /* parse latitude and longitude */
+ done |= pgl_scan(&strptr, &lat, &lon);
+ done |= pgl_scan(&strptr, &lat, &lon);
+ /* parse radius (while incrementing strptr by number of bytes parsed) */
+ valid = 0;
+ if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
+ }
+ /* require end of string and both latitude and longitude being parsed */
+ if (strptr[0] || done != PGL_SCAN_LATLON) {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+ errmsg("invalid input syntax for type ecircle: \"%s\"", str)
+ ));
+ }
+ /* allocate memory for result */
+ circle = (pgl_circle *)palloc(sizeof(pgl_circle));
+ /* set latitude, longitude, radius (and perform checks) */
+ pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
+ /* return result */
+ PG_RETURN_POINTER(circle);
+}
+
+/* parse cluster ("ecluster" in SQL) */
+PG_FUNCTION_INFO_V1(pgl_ecluster_in);
+Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
+ int i;
+ char *str = PG_GETARG_CSTRING(0); /* input string */
+ char *str_lower; /* lower case version of input string */
+ char *strptr; /* pointer to current reading position of input */
+ int npoints_total = 0; /* total number of points in cluster */
+ int nentries = 0; /* total number of entries */
+ pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */
+ int entries_buflen = 4; /* maximum number of elements in entries array */
+ int valid; /* number of valid chars processed */
+ double lat, lon; /* latitude and longitude of parsed point */
+ int entrytype; /* current entry type */
+ int npoints; /* number of points in current entry */
+ pgl_point *points; /* array of pgl_point for pgl_newentry */
+ int points_buflen; /* maximum number of elements in points array */
+ int done; /* return value of pgl_scan function */
+ pgl_cluster *cluster; /* created cluster */
+ /* lowercase input */
+ str_lower = psprintf("%s", str);
+ for (strptr=str_lower; *strptr; strptr++) {
+ if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
+ }
+ /* reset reading position to start of (lowercase) string */
+ strptr = str_lower;
+ /* allocate initial buffer for entries */
+ entries = palloc(entries_buflen * sizeof(pgl_newentry));
+ /* parse until end of string */
+ while (strptr[0]) {
+ /* require previous white-space or closing parenthesis before next token */
+ if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
+ goto pgl_ecluster_in_error;
+ }
+ /* ignore token "empty" */
+ valid = 0; sscanf(strptr, " empty %n", &valid);
+ if (valid) { strptr += valid; continue; }
+ /* test for "point" token */
+ valid = 0; sscanf(strptr, " point ( %n", &valid);
+ if (valid) {
+ strptr += valid;
+ entrytype = PGL_ENTRY_POINT;
+ goto pgl_ecluster_in_type_ok;
+ }
+ /* test for "path" token */
+ valid = 0; sscanf(strptr, " path ( %n", &valid);
+ if (valid) {
+ strptr += valid;
+ entrytype = PGL_ENTRY_PATH;
+ goto pgl_ecluster_in_type_ok;
+ }
+ /* test for "outline" token */
+ valid = 0; sscanf(strptr, " outline ( %n", &valid);
+ if (valid) {
+ strptr += valid;
+ entrytype = PGL_ENTRY_OUTLINE;
+ goto pgl_ecluster_in_type_ok;
+ }
+ /* test for "polygon" token */
+ valid = 0; sscanf(strptr, " polygon ( %n", &valid);
+ if (valid) {
+ strptr += valid;
+ entrytype = PGL_ENTRY_POLYGON;
+ goto pgl_ecluster_in_type_ok;
+ }
+ /* error if no valid token found */
+ goto pgl_ecluster_in_error;
+ pgl_ecluster_in_type_ok:
+ /* check if pgl_newentry array needs to grow */
+ if (nentries == entries_buflen) {
+ pgl_newentry *newbuf;
+ entries_buflen *= 2;
+ newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
+ memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
+ pfree(entries);
+ entries = newbuf;
+ }
+ /* reset number of points for current entry */
+ npoints = 0;
+ /* allocate array for points */
+ points_buflen = 4;
+ points = palloc(points_buflen * sizeof(pgl_point));
+ /* parse until closing parenthesis */
+ while (strptr[0] != ')') {
+ /* error on unexpected end of string */
+ if (strptr[0] == 0) goto pgl_ecluster_in_error;
+ /* mark neither latitude nor longitude as read */
+ done = PGL_SCAN_NONE;
+ /* require white-space before second, third, etc. point */
+ if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
+ /* scan latitude (or longitude) */
+ done |= pgl_scan(&strptr, &lat, &lon);
+ /* require white-space before second coordinate */
+ if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
+ /* scan longitude (or latitude) */
+ done |= pgl_scan(&strptr, &lat, &lon);
+ /* error unless both latitude and longitude were parsed */
+ if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
+ /* throw error if number of points is too high */
+ if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+ errmsg(
+ "too many points for ecluster entry (maximum %i)",
+ PGL_CLUSTER_MAXPOINTS
+ )
+ ));
+ }
+ /* check if pgl_point array needs to grow */
+ if (npoints == points_buflen) {
+ pgl_point *newbuf;
+ points_buflen *= 2;
+ newbuf = palloc(points_buflen * sizeof(pgl_point));
+ memcpy(newbuf, points, npoints * sizeof(pgl_point));
+ pfree(points);
+ points = newbuf;
+ }
+ /* append point to pgl_point array (includes checks) */
+ pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
+ /* increase total number of points */
+ npoints_total++;
+ }
+ /* error if entry has no points */
+ if (!npoints) goto pgl_ecluster_in_error;
+ /* entries with one point are automatically of type "point" */
+ if (npoints == 1) entrytype = PGL_ENTRY_POINT;
+ /* if entries have more than one point */
+ else {
+ /* throw error if entry type is "point" */
+ if (entrytype == PGL_ENTRY_POINT) {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+ errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
+ ));
+ }
+ /* coerce outlines and polygons with more than 2 points to be a path */
+ if (npoints == 2) entrytype = PGL_ENTRY_PATH;
+ }
+ /* append entry to pgl_newentry array */
+ entries[nentries].entrytype = entrytype;
+ entries[nentries].npoints = npoints;
+ entries[nentries].points = points;
+ nentries++;
+ /* consume closing parenthesis */
+ strptr++;
+ /* consume white-space */
+ while (isspace(strptr[0])) strptr++;
+ }
+ /* free lower case string */
+ pfree(str_lower);
+ /* create cluster from pgl_newentry array */
+ cluster = pgl_new_cluster(nentries, entries);
+ /* free pgl_newentry array */
+ for (i=0; ilat);
+ pgl_print_lon(lonstr, point->lon);
+ PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
+}
+
+/* convert box ("ebox") to string representation */
+PG_FUNCTION_INFO_V1(pgl_ebox_out);
+Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
+ double lon_min = box->lon_min;
+ double lon_max = box->lon_max;
+ char lat_min_str[PGL_NUMBUFLEN];
+ char lat_max_str[PGL_NUMBUFLEN];
+ char lon_min_str[PGL_NUMBUFLEN];
+ char lon_max_str[PGL_NUMBUFLEN];
+ /* return string "empty" if box is set to be empty */
+ if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
+ /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
+ /* (required since pgl_box_in orders the longitude boundaries) */
+ if (lon_min > lon_max) {
+ if (lon_min + lon_max >= 0) lon_min -= 360;
+ else lon_max += 360;
+ }
+ /* format and return result */
+ pgl_print_lat(lat_min_str, box->lat_min);
+ pgl_print_lat(lat_max_str, box->lat_max);
+ pgl_print_lon(lon_min_str, lon_min);
+ pgl_print_lon(lon_max_str, lon_max);
+ PG_RETURN_CSTRING(psprintf(
+ "%s %s %s %s",
+ lat_min_str, lon_min_str, lat_max_str, lon_max_str
+ ));
+}
+
+/* convert circle ("ecircle") to string representation */
+PG_FUNCTION_INFO_V1(pgl_ecircle_out);
+Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
+ char latstr[PGL_NUMBUFLEN];
+ char lonstr[PGL_NUMBUFLEN];
+ char radstr[PGL_NUMBUFLEN];
+ pgl_print_lat(latstr, circle->center.lat);
+ pgl_print_lon(lonstr, circle->center.lon);
+ pgl_print_float(radstr, circle->radius);
+ PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
+}
+
+/* convert cluster ("ecluster") to string representation */
+PG_FUNCTION_INFO_V1(pgl_ecluster_out);
+Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */
+ char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */
+ char ***strings; /* array of array of strings */
+ char *string; /* string of current token */
+ char *res, *resptr; /* result and pointer to current write position */
+ size_t reslen = 1; /* length of result (init with 1 for terminator) */
+ int npoints; /* number of points of current entry */
+ int i, j; /* i: entry, j: point in entry */
+ /* handle empty clusters */
+ if (cluster->nentries == 0) {
+ /* free detoasted cluster (if copy) */
+ PG_FREE_IF_COPY(cluster, 0);
+ /* return static result */
+ PG_RETURN_CSTRING("empty");
+ }
+ /* allocate array of array of strings */
+ strings = palloc(cluster->nentries * sizeof(char **));
+ /* iterate over all entries in cluster */
+ for (i=0; inentries; i++) {
+ /* get number of points in entry */
+ npoints = cluster->entries[i].npoints;
+ /* allocate array of strings (one string for each point plus two extra) */
+ strings[i] = palloc((2 + npoints) * sizeof(char *));
+ /* determine opening string */
+ switch (cluster->entries[i].entrytype) {
+ case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break;
+ case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break;
+ case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
+ case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
+ default: string = (i==0)?"unknown" :" unknown";
+ }
+ /* use opening string as first string in array */
+ strings[i][0] = string;
+ /* update result length (for allocating result string later) */
+ reslen += strlen(string);
+ /* iterate over all points */
+ for (j=0; jnentries; i++) {
+ npoints = cluster->entries[i].npoints;
+ for (j=0; jlat = pq_getmsgfloat8(buf);
+ point->lon = pq_getmsgfloat8(buf);
+ PG_RETURN_POINTER(point);
+}
+
+/* binary input function for box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_recv);
+Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
+ StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
+ pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
+ box->lat_min = pq_getmsgfloat8(buf);
+ box->lat_max = pq_getmsgfloat8(buf);
+ box->lon_min = pq_getmsgfloat8(buf);
+ box->lon_max = pq_getmsgfloat8(buf);
+ PG_RETURN_POINTER(box);
+}
+
+/* binary input function for circle ("ecircle") */
+PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
+Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
+ StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
+ pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
+ circle->center.lat = pq_getmsgfloat8(buf);
+ circle->center.lon = pq_getmsgfloat8(buf);
+ circle->radius = pq_getmsgfloat8(buf);
+ PG_RETURN_POINTER(circle);
+}
+
+/* TODO: binary receive function for cluster */
+
+/* binary output function for point ("epoint") */
+PG_FUNCTION_INFO_V1(pgl_epoint_send);
+Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ StringInfoData buf;
+ pq_begintypsend(&buf);
+ pq_sendfloat8(&buf, point->lat);
+ pq_sendfloat8(&buf, point->lon);
+ PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
+}
+
+/* binary output function for box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_send);
+Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
+ StringInfoData buf;
+ pq_begintypsend(&buf);
+ pq_sendfloat8(&buf, box->lat_min);
+ pq_sendfloat8(&buf, box->lat_max);
+ pq_sendfloat8(&buf, box->lon_min);
+ pq_sendfloat8(&buf, box->lon_max);
+ PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
+}
+
+/* binary output function for circle ("ecircle") */
+PG_FUNCTION_INFO_V1(pgl_ecircle_send);
+Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
+ StringInfoData buf;
+ pq_begintypsend(&buf);
+ pq_sendfloat8(&buf, circle->center.lat);
+ pq_sendfloat8(&buf, circle->center.lon);
+ pq_sendfloat8(&buf, circle->radius);
+ PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
+}
+
+/* TODO: binary send functions for cluster */
+
+/* cast point ("epoint") to box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
+Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_box *box = palloc(sizeof(pgl_box));
+ box->lat_min = point->lat;
+ box->lat_max = point->lat;
+ box->lon_min = point->lon;
+ box->lon_max = point->lon;
+ PG_RETURN_POINTER(box);
+}
+
+/* cast point ("epoint") to circle ("ecircle") */
+PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
+Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_circle *circle = palloc(sizeof(pgl_box));
+ circle->center = *point;
+ circle->radius = 0;
+ PG_RETURN_POINTER(circle);
+}
+
+/* cast point ("epoint") to cluster ("ecluster") */
+PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
+Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_newentry entry;
+ pgl_cluster *cluster;
+ entry.entrytype = PGL_ENTRY_POINT;
+ entry.npoints = 1;
+ entry.points = point;
+ cluster = pgl_new_cluster(1, &entry);
+ pgl_finalize_cluster(cluster); /* NOTE: should not fail */
+ PG_RETURN_POINTER(cluster);
+}
+
+/* cast box ("ebox") to cluster ("ecluster") */
+#define pgl_ebox_to_ecluster_macro(i, a, b) \
+ entries[i].entrytype = PGL_ENTRY_POLYGON; \
+ entries[i].npoints = 4; \
+ entries[i].points = points[i]; \
+ points[i][0].lat = box->lat_min; \
+ points[i][0].lon = (a); \
+ points[i][1].lat = box->lat_min; \
+ points[i][1].lon = (b); \
+ points[i][2].lat = box->lat_max; \
+ points[i][2].lon = (b); \
+ points[i][3].lat = box->lat_max; \
+ points[i][3].lon = (a);
+PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
+Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
+ double lon, dlon;
+ int nentries;
+ pgl_newentry entries[3];
+ pgl_point points[3][4];
+ pgl_cluster *cluster;
+ if (box->lat_min > box->lat_max) {
+ nentries = 0;
+ } else if (box->lon_min > box->lon_max) {
+ if (box->lon_min < 0) {
+ lon = pgl_round((box->lon_min + 180) / 2.0);
+ nentries = 3;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
+ pgl_ebox_to_ecluster_macro(1, lon, 180);
+ pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
+ } else if (box->lon_max > 0) {
+ lon = pgl_round((box->lon_max - 180) / 2.0);
+ nentries = 3;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
+ pgl_ebox_to_ecluster_macro(1, -180, lon);
+ pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
+ } else {
+ nentries = 2;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
+ pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
+ }
+ } else {
+ dlon = pgl_round(box->lon_max - box->lon_min);
+ if (dlon < 180) {
+ nentries = 1;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
+ } else {
+ lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
+ if (
+ pgl_round(lon - box->lon_min) < 180 &&
+ pgl_round(box->lon_max - lon) < 180
+ ) {
+ nentries = 2;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
+ pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
+ } else {
+ nentries = 3;
+ pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
+ pgl_ebox_to_ecluster_macro(1, -60, 60);
+ pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
+ }
+ }
+ }
+ cluster = pgl_new_cluster(nentries, entries);
+ pgl_finalize_cluster(cluster); /* NOTE: should not fail */
+ PG_RETURN_POINTER(cluster);
+}
+
+/* extract latitude from point ("epoint") */
+PG_FUNCTION_INFO_V1(pgl_epoint_lat);
+Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
+}
+
+/* extract longitude from point ("epoint") */
+PG_FUNCTION_INFO_V1(pgl_epoint_lon);
+Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
+}
+
+/* extract minimum latitude from box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
+Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
+}
+
+/* extract maximum latitude from box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
+Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
+}
+
+/* extract minimum longitude from box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
+Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
+}
+
+/* extract maximum longitude from box ("ebox") */
+PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
+Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
+}
+
+/* extract center point from circle ("ecircle") */
+PG_FUNCTION_INFO_V1(pgl_ecircle_center);
+Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
+ PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
+}
+
+/* extract radius from circle ("ecircle") */
+PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
+Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
+ PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
+}
+
+/* check if point is inside box (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
+Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
+ PG_RETURN_BOOL(pgl_point_in_box(point, box));
+}
+
+/* check if point is inside circle (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
+Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
+ PG_RETURN_BOOL(
+ pgl_distance(
+ point->lat, point->lon,
+ circle->center.lat, circle->center.lon
+ ) <= circle->radius
+ );
+}
+
+/* check if point is inside cluster (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
+Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval;
+ /* points outside bounding circle are always assumed to be non-overlapping
+ (necessary for consistent table and index scans) */
+ if (
+ pgl_distance(
+ point->lat, point->lon,
+ cluster->bounding.center.lat, cluster->bounding.center.lon
+ ) > cluster->bounding.radius
+ ) retval = false;
+ else retval = pgl_point_in_cluster(point, cluster, false);
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
+Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval = pgl_distance(
+ point->lat, point->lon,
+ cluster->bounding.center.lat, cluster->bounding.center.lon
+ ) <= cluster->bounding.radius;
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if two boxes overlap (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
+Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
+ pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
+ pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
+ PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
+}
+
+/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
+Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
+ PG_RETURN_BOOL(
+ pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
+ );
+}
+
+/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
+Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
+ pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval = pgl_estimate_point_box_distance(
+ &cluster->bounding.center,
+ box
+ ) <= cluster->bounding.radius;
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if two circles overlap (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
+Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
+ pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
+ pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
+ PG_RETURN_BOOL(
+ pgl_distance(
+ circle1->center.lat, circle1->center.lon,
+ circle2->center.lat, circle2->center.lon
+ ) <= circle1->radius + circle2->radius
+ );
+}
+
+/* check if circle and cluster overlap (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
+Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval = (
+ pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
+ );
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
+Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval = pgl_distance(
+ circle->center.lat, circle->center.lon,
+ cluster->bounding.center.lat, cluster->bounding.center.lon
+ ) <= circle->radius + cluster->bounding.radius;
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if two clusters overlap (overlap operator "&&") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
+Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval;
+ /* clusters with non-touching bounding circles are always assumed to be
+ non-overlapping (improves performance and is necessary for consistent
+ table and index scans) */
+ if (
+ pgl_distance(
+ cluster1->bounding.center.lat, cluster1->bounding.center.lon,
+ cluster2->bounding.center.lat, cluster2->bounding.center.lon
+ ) > cluster1->bounding.radius + cluster2->bounding.radius
+ ) retval = false;
+ else retval = pgl_clusters_overlap(cluster1, cluster2);
+ PG_FREE_IF_COPY(cluster1, 0);
+ PG_FREE_IF_COPY(cluster2, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
+Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval = pgl_distance(
+ cluster1->bounding.center.lat, cluster1->bounding.center.lon,
+ cluster2->bounding.center.lat, cluster2->bounding.center.lon
+ ) <= cluster1->bounding.radius + cluster2->bounding.radius;
+ PG_FREE_IF_COPY(cluster1, 0);
+ PG_FREE_IF_COPY(cluster2, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* check if second cluster is in first cluster (cont. operator "@>) in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
+Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
+ pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ bool retval;
+ /* clusters with non-touching bounding circles are always assumed to be
+ non-overlapping (improves performance and is necessary for consistent
+ table and index scans) */
+ if (
+ pgl_distance(
+ outer->bounding.center.lat, outer->bounding.center.lon,
+ inner->bounding.center.lat, inner->bounding.center.lon
+ ) > outer->bounding.radius + inner->bounding.radius
+ ) retval = false;
+ else retval = pgl_cluster_in_cluster(outer, inner);
+ PG_FREE_IF_COPY(outer, 0);
+ PG_FREE_IF_COPY(inner, 1);
+ PG_RETURN_BOOL(retval);
+}
+
+/* calculate distance between two points ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_distance);
+Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
+ pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
+ PG_RETURN_FLOAT8(pgl_distance(
+ point1->lat, point1->lon, point2->lat, point2->lon
+ ));
+}
+
+/* calculate point to circle distance ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
+Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
+ double distance = pgl_distance(
+ point->lat, point->lon, circle->center.lat, circle->center.lon
+ ) - circle->radius;
+ PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
+}
+
+/* calculate point to cluster distance ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
+Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ double distance = pgl_point_cluster_distance(point, cluster);
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_FLOAT8(distance);
+}
+
+/* calculate distance between two circles ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
+Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
+ pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
+ pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
+ double distance = pgl_distance(
+ circle1->center.lat, circle1->center.lon,
+ circle2->center.lat, circle2->center.lon
+ ) - (circle1->radius + circle2->radius);
+ PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
+}
+
+/* calculate circle to cluster distance ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
+Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
+ pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ double distance = (
+ pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
+ );
+ PG_FREE_IF_COPY(cluster, 1);
+ PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
+}
+
+/* calculate distance between two clusters ("<->" operator) in SQL */
+PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
+Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ double retval = pgl_cluster_distance(cluster1, cluster2);
+ PG_FREE_IF_COPY(cluster1, 0);
+ PG_FREE_IF_COPY(cluster2, 1);
+ PG_RETURN_FLOAT8(retval);
+}
+
+PG_FUNCTION_INFO_V1(pgl_ecluster_monte_carlo_area);
+Datum pgl_ecluster_monte_carlo_area(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ int samples = PG_GETARG_INT32(1);
+ double retval = pgl_monte_carlo_area(cluster, samples);
+ PG_FREE_IF_COPY(cluster, 0);
+ PG_RETURN_FLOAT8(retval);
+}
+
+PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_fair_distance);
+Datum pgl_ecluster_epoint_fair_distance(PG_FUNCTION_ARGS) {
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ pgl_point *point = (pgl_point *)PG_GETARG_POINTER(1);
+ int samples = PG_GETARG_INT32(2);
+ double retval = pgl_fair_distance(point, cluster, samples);
+ PG_FREE_IF_COPY(cluster, 0);
+ PG_RETURN_FLOAT8(retval);
+}
+
+
+/*-----------------------------------------------------------*
+ * B-tree comparison operators and index support functions *
+ *-----------------------------------------------------------*/
+
+/* macro for a B-tree operator (without detoasting) */
+#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
+ PG_FUNCTION_INFO_V1(func); \
+ Datum func(PG_FUNCTION_ARGS) { \
+ type *a = (type *)PG_GETARG_POINTER(0); \
+ type *b = (type *)PG_GETARG_POINTER(1); \
+ PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
+ }
+
+/* macro for a B-tree comparison function (without detoasting) */
+#define PGL_BTREE_CMP(func, type, cmpfunc) \
+ PG_FUNCTION_INFO_V1(func); \
+ Datum func(PG_FUNCTION_ARGS) { \
+ type *a = (type *)PG_GETARG_POINTER(0); \
+ type *b = (type *)PG_GETARG_POINTER(1); \
+ PG_RETURN_INT32(cmpfunc(a, b)); \
+ }
+
+/* macro for a B-tree operator (with detoasting) */
+#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
+ PG_FUNCTION_INFO_V1(func); \
+ Datum func(PG_FUNCTION_ARGS) { \
+ bool res; \
+ type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
+ type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
+ res = cmpfunc(a, b) oper 0; \
+ PG_FREE_IF_COPY(a, 0); \
+ PG_FREE_IF_COPY(b, 1); \
+ PG_RETURN_BOOL(res); \
+ }
+
+/* macro for a B-tree comparison function (with detoasting) */
+#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
+ PG_FUNCTION_INFO_V1(func); \
+ Datum func(PG_FUNCTION_ARGS) { \
+ int32_t res; \
+ type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
+ type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
+ res = cmpfunc(a, b); \
+ PG_FREE_IF_COPY(a, 0); \
+ PG_FREE_IF_COPY(b, 1); \
+ PG_RETURN_INT32(res); \
+ }
+
+/* B-tree operators and comparison function for point */
+PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
+PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
+PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
+PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
+PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
+PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
+PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
+
+/* B-tree operators and comparison function for box */
+PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
+PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
+PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
+PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
+PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
+PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
+PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
+
+/* B-tree operators and comparison function for circle */
+PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
+PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
+PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
+PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
+PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
+PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
+PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
+
+
+/*--------------------------------*
+ * GiST index support functions *
+ *--------------------------------*/
+
+/* GiST "consistent" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_consistent);
+Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
+ GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+ pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
+ StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
+ bool *recheck = (bool *)PG_GETARG_POINTER(4);
+ /* demand recheck because index and query methods are lossy */
+ *recheck = true;
+ /* strategy number aliases for different operators using the same strategy */
+ strategy %= 100;
+ /* strategy number 11: equality of two points */
+ if (strategy == 11) {
+ /* query datum is another point */
+ pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
+ /* convert other point to key */
+ pgl_pointkey querykey;
+ pgl_point_to_key(query, querykey);
+ /* return true if both keys overlap */
+ PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
+ }
+ /* strategy number 13: equality of two circles */
+ if (strategy == 13) {
+ /* query datum is another circle */
+ pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
+ /* convert other circle to key */
+ pgl_areakey querykey;
+ pgl_circle_to_key(query, querykey);
+ /* return true if both keys overlap */
+ PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
+ }
+ /* for all remaining strategies, keys on empty objects produce no match */
+ /* (check necessary because query radius may be infinite) */
+ if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
+ /* strategy number 21: overlapping with point */
+ if (strategy == 21) {
+ /* query datum is a point */
+ pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
+ /* return true if estimated distance (allowed to be smaller than real
+ distance) between index key and point is zero */
+ PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
+ }
+ /* strategy number 22: (point) overlapping with box */
+ if (strategy == 22) {
+ /* query datum is a box */
+ pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
+ /* determine bounding box of indexed key */
+ pgl_box keybox;
+ pgl_key_to_box(key, &keybox);
+ /* return true if query box overlaps with bounding box of indexed key */
+ PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
+ }
+ /* strategy number 23: overlapping with circle */
+ if (strategy == 23) {
+ /* query datum is a circle */
+ pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
+ /* return true if estimated distance (allowed to be smaller than real
+ distance) between index key and circle center is smaller than radius */
+ PG_RETURN_BOOL(
+ pgl_estimate_key_distance(key, &(query->center)) <= query->radius
+ );
+ }
+ /* strategy number 24: overlapping with cluster */
+ if (strategy == 24) {
+ bool retval; /* return value */
+ /* query datum is a cluster */
+ pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ /* return true if estimated distance (allowed to be smaller than real
+ distance) between index key and circle center is smaller than radius */
+ retval = (
+ pgl_estimate_key_distance(key, &(query->bounding.center)) <=
+ query->bounding.radius
+ );
+ PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */
+ PG_RETURN_BOOL(retval);
+ }
+ /* throw error for any unknown strategy number */
+ elog(ERROR, "unrecognized strategy number: %d", strategy);
+}
+
+/* GiST "union" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_union);
+Datum pgl_gist_union(PG_FUNCTION_ARGS) {
+ GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
+ pgl_keyptr out; /* return value (to be palloc'ed) */
+ int i;
+ /* determine key size */
+ size_t keysize = PGL_KEY_IS_AREAKEY(
+ (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
+ ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
+ /* begin with first key as result */
+ out = palloc(keysize);
+ memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
+ /* unite current result with second, third, etc. key */
+ for (i=1; in; i++) {
+ pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
+ }
+ /* return result */
+ PG_RETURN_POINTER(out);
+}
+
+/* GiST "compress" support function for indicis on points */
+PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
+Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
+ GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+ GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
+ /* only transform new leaves */
+ if (entry->leafkey) {
+ /* get point to be transformed */
+ pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
+ /* allocate memory for key */
+ pgl_keyptr key = palloc(sizeof(pgl_pointkey));
+ /* transform point to key */
+ pgl_point_to_key(point, key);
+ /* create new GISTENTRY structure as return value */
+ retval = palloc(sizeof(GISTENTRY));
+ gistentryinit(
+ *retval, PointerGetDatum(key),
+ entry->rel, entry->page, entry->offset, FALSE
+ );
+ } else {
+ /* inner nodes have already been transformed */
+ retval = entry;
+ }
+ /* return pointer to old or new GISTENTRY structure */
+ PG_RETURN_POINTER(retval);
+}
+
+/* GiST "compress" support function for indicis on circles */
+PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
+Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
+ GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+ GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
+ /* only transform new leaves */
+ if (entry->leafkey) {
+ /* get circle to be transformed */
+ pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
+ /* allocate memory for key */
+ pgl_keyptr key = palloc(sizeof(pgl_areakey));
+ /* transform circle to key */
+ pgl_circle_to_key(circle, key);
+ /* create new GISTENTRY structure as return value */
+ retval = palloc(sizeof(GISTENTRY));
+ gistentryinit(
+ *retval, PointerGetDatum(key),
+ entry->rel, entry->page, entry->offset, FALSE
+ );
+ } else {
+ /* inner nodes have already been transformed */
+ retval = entry;
+ }
+ /* return pointer to old or new GISTENTRY structure */
+ PG_RETURN_POINTER(retval);
+}
+
+/* GiST "compress" support function for indices on clusters */
+PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
+Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
+ GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+ GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */
+ /* only transform new leaves */
+ if (entry->leafkey) {
+ /* get cluster to be transformed (detoasting necessary!) */
+ pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
+ /* allocate memory for key */
+ pgl_keyptr key = palloc(sizeof(pgl_areakey));
+ /* transform cluster to key */
+ pgl_circle_to_key(&(cluster->bounding), key);
+ /* create new GISTENTRY structure as return value */
+ retval = palloc(sizeof(GISTENTRY));
+ gistentryinit(
+ *retval, PointerGetDatum(key),
+ entry->rel, entry->page, entry->offset, FALSE
+ );
+ /* free detoasted datum */
+ if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
+ } else {
+ /* inner nodes have already been transformed */
+ retval = entry;
+ }
+ /* return pointer to old or new GISTENTRY structure */
+ PG_RETURN_POINTER(retval);
+}
+
+/* GiST "decompress" support function for indices */
+PG_FUNCTION_INFO_V1(pgl_gist_decompress);
+Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
+ /* return passed pointer without transformation */
+ PG_RETURN_POINTER(PG_GETARG_POINTER(0));
+}
+
+/* GiST "penalty" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_penalty);
+Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
+ GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
+ GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
+ float *penalty = (float *)PG_GETARG_POINTER(2);
+ /* get original key and key to insert */
+ pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
+ pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
+ /* copy original key */
+ union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
+ if (PGL_KEY_IS_AREAKEY(orig)) {
+ memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
+ } else {
+ memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
+ }
+ /* calculate union of both keys */
+ pgl_unite_keys((pgl_keyptr)&union_key, new);
+ /* penalty equal to reduction of key length (logarithm of added area) */
+ /* (return value by setting referenced value and returning pointer) */
+ *penalty = (
+ PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
+ );
+ PG_RETURN_POINTER(penalty);
+}
+
+/* GiST "picksplit" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
+Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
+ GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
+ GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
+ OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */
+ union {
+ pgl_pointkey pointkey;
+ pgl_areakey areakey;
+ } union_all; /* union of all keys (to be calculated from scratch)
+ (later cut in half) */
+ int is_areakey = PGL_KEY_IS_AREAKEY(
+ (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
+ );
+ int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
+ pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */
+ pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */
+ pgl_keyptr key; /* current key to be processed */
+ /* allocate memory for array of left and right keys, set counts to zero */
+ v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
+ v->spl_nleft = 0;
+ v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
+ v->spl_nright = 0;
+ /* calculate union of all keys from scratch */
+ memcpy(
+ (pgl_keyptr)&union_all,
+ (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
+ keysize
+ );
+ for (i=FirstOffsetNumber+1; in; i=OffsetNumberNext(i)) {
+ pgl_unite_keys(
+ (pgl_keyptr)&union_all,
+ (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
+ );
+ }
+ /* check if trivial split is necessary due to exhausted key length */
+ /* (Note: keys for empty objects must have node depth set to maximum) */
+ if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
+ is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
+ )) {
+ /* half of all keys go left */
+ for (
+ i=FirstOffsetNumber;
+ in - FirstOffsetNumber)/2;
+ i=OffsetNumberNext(i)
+ ) {
+ /* pointer to current key */
+ key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
+ /* update unionL */
+ /* check if key is first key that goes left */
+ if (!v->spl_nleft) {
+ /* first key that goes left is just copied to unionL */
+ memcpy(unionL, key, keysize);
+ } else {
+ /* unite current value and next key */
+ pgl_unite_keys(unionL, key);
+ }
+ /* append offset number to list of keys that go left */
+ v->spl_left[v->spl_nleft++] = i;
+ }
+ /* other half goes right */
+ for (
+ i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
+ in;
+ i=OffsetNumberNext(i)
+ ) {
+ /* pointer to current key */
+ key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
+ /* update unionR */
+ /* check if key is first key that goes right */
+ if (!v->spl_nright) {
+ /* first key that goes right is just copied to unionR */
+ memcpy(unionR, key, keysize);
+ } else {
+ /* unite current value and next key */
+ pgl_unite_keys(unionR, key);
+ }
+ /* append offset number to list of keys that go right */
+ v->spl_right[v->spl_nright++] = i;
+ }
+ }
+ /* otherwise, a non-trivial split is possible */
+ else {
+ /* cut covered area in half */
+ /* (union_all then refers to area of keys that go left) */
+ /* check if union of all keys covers empty and non-empty objects */
+ if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
+ /* if yes, split into empty and non-empty objects */
+ pgl_key_set_empty((pgl_keyptr)&union_all);
+ } else {
+ /* otherwise split by next bit */
+ ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
+ /* NOTE: type bit conserved */
+ }
+ /* determine for each key if it goes left or right */
+ for (i=FirstOffsetNumber; in; i=OffsetNumberNext(i)) {
+ /* pointer to current key */
+ key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
+ /* keys within one half of the area go left */
+ if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
+ /* update unionL */
+ /* check if key is first key that goes left */
+ if (!v->spl_nleft) {
+ /* first key that goes left is just copied to unionL */
+ memcpy(unionL, key, keysize);
+ } else {
+ /* unite current value of unionL and processed key */
+ pgl_unite_keys(unionL, key);
+ }
+ /* append offset number to list of keys that go left */
+ v->spl_left[v->spl_nleft++] = i;
+ }
+ /* the other keys go right */
+ else {
+ /* update unionR */
+ /* check if key is first key that goes right */
+ if (!v->spl_nright) {
+ /* first key that goes right is just copied to unionR */
+ memcpy(unionR, key, keysize);
+ } else {
+ /* unite current value of unionR and processed key */
+ pgl_unite_keys(unionR, key);
+ }
+ /* append offset number to list of keys that go right */
+ v->spl_right[v->spl_nright++] = i;
+ }
+ }
+ }
+ /* store unions in return value */
+ v->spl_ldatum = PointerGetDatum(unionL);
+ v->spl_rdatum = PointerGetDatum(unionR);
+ /* return all results */
+ PG_RETURN_POINTER(v);
+}
+
+/* GiST "same"/"equal" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_same);
+Datum pgl_gist_same(PG_FUNCTION_ARGS) {
+ pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
+ pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
+ bool *boolptr = (bool *)PG_GETARG_POINTER(2);
+ /* two keys are equal if they are binary equal */
+ /* (return result by setting referenced boolean and returning pointer) */
+ *boolptr = !memcmp(
+ key1,
+ key2,
+ PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
+ );
+ PG_RETURN_POINTER(boolptr);
+}
+
+/* GiST "distance" support function */
+PG_FUNCTION_INFO_V1(pgl_gist_distance);
+Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
+ GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
+ pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
+ StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
+ bool *recheck = (bool *)PG_GETARG_POINTER(4);
+ double distance; /* return value */
+ /* demand recheck because distance is just an estimation */
+ /* (real distance may be bigger) */
+ *recheck = true;
+ /* strategy number aliases for different operators using the same strategy */
+ strategy %= 100;
+ /* strategy number 31: distance to point */
+ if (strategy == 31) {
+ /* query datum is a point */
+ pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
+ /* use pgl_estimate_pointkey_distance() function to compute result */
+ distance = pgl_estimate_key_distance(key, query);
+ /* avoid infinity (reserved!) */
+ if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
+ /* return result */
+ PG_RETURN_FLOAT8(distance);
+ }
+ /* strategy number 33: distance to circle */
+ if (strategy == 33) {
+ /* query datum is a circle */
+ pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
+ /* estimate distance to circle center and substract circle radius */
+ distance = (
+ pgl_estimate_key_distance(key, &(query->center)) - query->radius
+ );
+ /* convert non-positive values to zero and avoid infinity (reserved!) */
+ if (distance <= 0) distance = 0;
+ else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
+ /* return result */
+ PG_RETURN_FLOAT8(distance);
+ }
+ /* strategy number 34: distance to cluster */
+ if (strategy == 34) {
+ /* query datum is a cluster */
+ pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+ /* estimate distance to bounding center and substract bounding radius */
+ distance = (
+ pgl_estimate_key_distance(key, &(query->bounding.center)) -
+ query->bounding.radius
+ );
+ /* convert non-positive values to zero and avoid infinity (reserved!) */
+ if (distance <= 0) distance = 0;
+ else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
+ /* free detoasted cluster (if copy) */
+ PG_FREE_IF_COPY(query, 1);
+ /* return result */
+ PG_RETURN_FLOAT8(distance);
+ }
+ /* throw error for any unknown strategy number */
+ elog(ERROR, "unrecognized strategy number: %d", strategy);
+}
+
diff -r a90c0d750405 -r 1b9cd45e9e48 latlon.control
--- a/latlon.control Fri Oct 21 12:57:46 2016 +0200
+++ b/latlon.control Tue Oct 25 18:44:43 2016 +0200
@@ -1,3 +1,3 @@
comment = 'Geographic data types and spatial indexing for the WGS-84 spheroid'
-default_version = '0.8'
+default_version = '0.9'
relocatable = true