pgLatLon
changeset 42:1b9cd45e9e48
Bugfix for type casts to ecluster; New "fair_distance" function
author | jbe |
---|---|
date | Tue Oct 25 18:44:43 2016 +0200 (2016-10-25) |
parents | a90c0d750405 |
children | 3293906ea6a0 |
files | GNUmakefile README.html README.mkd latlon--0.9.sql latlon-v0008.c latlon.control |
line diff
1.1 --- a/GNUmakefile Fri Oct 21 12:57:46 2016 +0200 1.2 +++ b/GNUmakefile Tue Oct 25 18:44:43 2016 +0200 1.3 @@ -1,6 +1,6 @@ 1.4 EXTENSION = latlon 1.5 -DATA = latlon--0.7--0.8.sql latlon--0.8.sql 1.6 -MODULES = latlon-v0007 1.7 +DATA = latlon--0.7--0.8.sql latlon--0.8.sql latlon--0.9.sql 1.8 +MODULES = latlon-v0007 latlon-v0008 1.9 1.10 PG_CONFIG = pg_config 1.11 PGXS := $(shell $(PG_CONFIG) --pgxs)
2.1 --- a/README.html Fri Oct 21 12:57:46 2016 +0200 2.2 +++ b/README.html Tue Oct 25 18:44:43 2016 +0200 2.3 @@ -1,5 +1,5 @@ 2.4 -<html><head><title>pgLatLon v0.8 documentation</title></head><body> 2.5 -<h1>pgLatLon v0.8 documentation</h1> 2.6 +<html><head><title>pgLatLon v0.9 documentation</title></head><body> 2.7 +<h1>pgLatLon v0.9 documentation</h1> 2.8 2.9 <p>pgLatLon is a spatial database extension for the PostgreSQL object-relational 2.10 database management system providing geographic data types and spatial indexing 2.11 @@ -40,8 +40,8 @@ 2.12 <p>It is also possible to compile and install the extension without GNU Make as 2.13 follows:</p> 2.14 2.15 -<pre><code>cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c 2.16 -cp latlon-v0007.so `pg_config --pkglibdir` 2.17 +<pre><code>cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c 2.18 +cp latlon-v0008.so `pg_config --pkglibdir` 2.19 cp latlon.control `pg_config --sharedir`/extension/ 2.20 cp latlon--*.sql `pg_config --sharedir`/extension/ 2.21 </code></pre> 2.22 @@ -486,6 +486,40 @@ 2.23 2.24 <p>Same as <code>epoint(float8, float8)</code> but with arguments reversed.</p> 2.25 2.26 +<h4><code>fair_distance(ecluster, epoint,</code> samples <code>int4 = 10000)</code></h4> 2.27 + 2.28 +<p>When working with user-generated content, users may be tempted to create 2.29 +intentionally oversized objects in order to optimize search results in an 2.30 +unfair manner. The <code>fair_distance</code> function aims to handle this by returning an 2.31 +adjusted distance (i.e. distance increased by a penalty) if a geographic object 2.32 +(the <code>ecluster</code>) consists of more than one point.</p> 2.33 + 2.34 +<p>The first argument to this function is an <code>ecluster</code>, the second argument is a 2.35 +search point (<code>epoint</code>), and the third argument is an interger related to the 2.36 +precision (higher precision will require more computation time).</p> 2.37 + 2.38 +<p>The penalty by which the returned distance is increased fulfills (at least) the 2.39 +following properties:</p> 2.40 + 2.41 +<ul> 2.42 +<li>For search points far away from the <code>ecluster</code> (i.e. distance approaching 2.43 +infinity), the penalty approaches zero (i.e. <code>fair_distance</code> behaves the same 2.44 +as <code>distance</code>).</li> 2.45 +<li>If the <code>ecluster</code> consists of a set of points, the penalty for a search point 2.46 +close to one of that points (closer than half of the minimum distance between 2.47 +each pair of points in the <code>ecluster</code>) is chosen in such a way that the 2.48 +adjusted distance is equal to the distance from the search point to the 2.49 +closest point in the <code>ecluster</code> multiplied by the square root of the count of 2.50 +points in the <code>ecluster</code>.</li> 2.51 +</ul> 2.52 + 2.53 +<p>The function interally uses a Monte Carlo simulation to compute the result. The 2.54 +third parameter (which defaults to 10000) can be used to adjust the number of 2.55 +samples taken. It is ensured that the penalty is always positive, i.e. results 2.56 +returned by the <code>fair_distance</code> function are always equal to or greater than 2.57 +the results returned by the <code>distance</code> function regardless of stochastic 2.58 +effects.</p> 2.59 + 2.60 <h4><code>GeoJSON_to_epoint(jsonb, text)</code></h4> 2.61 2.62 <p>Maps a GeoJSON object of type "Point" or "Feature" (which contains a
3.1 --- a/README.mkd Fri Oct 21 12:57:46 2016 +0200 3.2 +++ b/README.mkd Tue Oct 25 18:44:43 2016 +0200 3.3 @@ -1,4 +1,4 @@ 3.4 -pgLatLon v0.8 documentation 3.5 +pgLatLon v0.9 documentation 3.6 =========================== 3.7 3.8 pgLatLon is a spatial database extension for the PostgreSQL object-relational 3.9 @@ -39,8 +39,8 @@ 3.10 It is also possible to compile and install the extension without GNU Make as 3.11 follows: 3.12 3.13 - cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c 3.14 - cp latlon-v0007.so `pg_config --pkglibdir` 3.15 + cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c 3.16 + cp latlon-v0008.so `pg_config --pkglibdir` 3.17 cp latlon.control `pg_config --sharedir`/extension/ 3.18 cp latlon--*.sql `pg_config --sharedir`/extension/ 3.19 3.20 @@ -469,6 +469,38 @@ 3.21 3.22 Same as `epoint(float8, float8)` but with arguments reversed. 3.23 3.24 +#### `fair_distance(ecluster, epoint,` samples `int4 = 10000)` 3.25 + 3.26 +When working with user-generated content, users may be tempted to create 3.27 +intentionally oversized objects in order to optimize search results in an 3.28 +unfair manner. The `fair_distance` function aims to handle this by returning an 3.29 +adjusted distance (i.e. distance increased by a penalty) if a geographic object 3.30 +(the `ecluster`) consists of more than one point. 3.31 + 3.32 +The first argument to this function is an `ecluster`, the second argument is a 3.33 +search point (`epoint`), and the third argument is an interger related to the 3.34 +precision (higher precision will require more computation time). 3.35 + 3.36 +The penalty by which the returned distance is increased fulfills (at least) the 3.37 +following properties: 3.38 + 3.39 +* For search points far away from the `ecluster` (i.e. distance approaching 3.40 + infinity), the penalty approaches zero (i.e. `fair_distance` behaves the same 3.41 + as `distance`). 3.42 +* If the `ecluster` consists of a set of points, the penalty for a search point 3.43 + close to one of that points (closer than half of the minimum distance between 3.44 + each pair of points in the `ecluster`) is chosen in such a way that the 3.45 + adjusted distance is equal to the distance from the search point to the 3.46 + closest point in the `ecluster` multiplied by the square root of the count of 3.47 + points in the `ecluster`. 3.48 + 3.49 +The function interally uses a Monte Carlo simulation to compute the result. The 3.50 +third parameter (which defaults to 10000) can be used to adjust the number of 3.51 +samples taken. It is ensured that the penalty is always positive, i.e. results 3.52 +returned by the `fair_distance` function are always equal to or greater than 3.53 +the results returned by the `distance` function regardless of stochastic 3.54 +effects. 3.55 + 3.56 #### `GeoJSON_to_epoint(jsonb, text)` 3.57 3.58 Maps a GeoJSON object of type "Point" or "Feature" (which contains a
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 4.2 +++ b/latlon--0.9.sql Tue Oct 25 18:44:43 2016 +0200 4.3 @@ -0,0 +1,1692 @@ 4.4 + 4.5 +---------------------------------------- 4.6 +-- forward declarations (shell types) -- 4.7 +---------------------------------------- 4.8 + 4.9 +CREATE TYPE epoint; 4.10 +CREATE TYPE ebox; 4.11 +CREATE TYPE ecircle; 4.12 +CREATE TYPE ecluster; 4.13 + 4.14 + 4.15 +------------------------------------------------------------ 4.16 +-- dummy input/output functions for dummy index key types -- 4.17 +------------------------------------------------------------ 4.18 + 4.19 +CREATE FUNCTION ekey_point_in_dummy(cstring) 4.20 + RETURNS ekey_point 4.21 + LANGUAGE C IMMUTABLE STRICT 4.22 + AS '$libdir/latlon-v0008', 'pgl_notimpl'; 4.23 + 4.24 +CREATE FUNCTION ekey_point_out_dummy(ekey_point) 4.25 + RETURNS cstring 4.26 + LANGUAGE C IMMUTABLE STRICT 4.27 + AS '$libdir/latlon-v0008', 'pgl_notimpl'; 4.28 + 4.29 +CREATE FUNCTION ekey_area_in_dummy(cstring) 4.30 + RETURNS ekey_area 4.31 + LANGUAGE C IMMUTABLE STRICT 4.32 + AS '$libdir/latlon-v0008', 'pgl_notimpl'; 4.33 + 4.34 +CREATE FUNCTION ekey_area_out_dummy(ekey_area) 4.35 + RETURNS cstring 4.36 + LANGUAGE C IMMUTABLE STRICT 4.37 + AS '$libdir/latlon-v0008', 'pgl_notimpl'; 4.38 + 4.39 + 4.40 +-------------------------- 4.41 +-- text input functions -- 4.42 +-------------------------- 4.43 + 4.44 +CREATE FUNCTION epoint_in(cstring) 4.45 + RETURNS epoint 4.46 + LANGUAGE C IMMUTABLE STRICT 4.47 + AS '$libdir/latlon-v0008', 'pgl_epoint_in'; 4.48 + 4.49 +CREATE FUNCTION ebox_in(cstring) 4.50 + RETURNS ebox 4.51 + LANGUAGE C IMMUTABLE STRICT 4.52 + AS '$libdir/latlon-v0008', 'pgl_ebox_in'; 4.53 + 4.54 +CREATE FUNCTION ecircle_in(cstring) 4.55 + RETURNS ecircle 4.56 + LANGUAGE C IMMUTABLE STRICT 4.57 + AS '$libdir/latlon-v0008', 'pgl_ecircle_in'; 4.58 + 4.59 +CREATE FUNCTION ecluster_in(cstring) 4.60 + RETURNS ecluster 4.61 + LANGUAGE C IMMUTABLE STRICT 4.62 + AS '$libdir/latlon-v0008', 'pgl_ecluster_in'; 4.63 + 4.64 + 4.65 +--------------------------- 4.66 +-- text output functions -- 4.67 +--------------------------- 4.68 + 4.69 +CREATE FUNCTION epoint_out(epoint) 4.70 + RETURNS cstring 4.71 + LANGUAGE C IMMUTABLE STRICT 4.72 + AS '$libdir/latlon-v0008', 'pgl_epoint_out'; 4.73 + 4.74 +CREATE FUNCTION ebox_out(ebox) 4.75 + RETURNS cstring 4.76 + LANGUAGE C IMMUTABLE STRICT 4.77 + AS '$libdir/latlon-v0008', 'pgl_ebox_out'; 4.78 + 4.79 +CREATE FUNCTION ecircle_out(ecircle) 4.80 + RETURNS cstring 4.81 + LANGUAGE C IMMUTABLE STRICT 4.82 + AS '$libdir/latlon-v0008', 'pgl_ecircle_out'; 4.83 + 4.84 +CREATE FUNCTION ecluster_out(ecluster) 4.85 + RETURNS cstring 4.86 + LANGUAGE C IMMUTABLE STRICT 4.87 + AS '$libdir/latlon-v0008', 'pgl_ecluster_out'; 4.88 + 4.89 + 4.90 +-------------------------- 4.91 +-- binary I/O functions -- 4.92 +-------------------------- 4.93 + 4.94 +CREATE FUNCTION epoint_recv(internal) 4.95 + RETURNS epoint 4.96 + LANGUAGE C IMMUTABLE STRICT 4.97 + AS '$libdir/latlon-v0008', 'pgl_epoint_recv'; 4.98 + 4.99 +CREATE FUNCTION ebox_recv(internal) 4.100 + RETURNS ebox 4.101 + LANGUAGE C IMMUTABLE STRICT 4.102 + AS '$libdir/latlon-v0008', 'pgl_ebox_recv'; 4.103 + 4.104 +CREATE FUNCTION ecircle_recv(internal) 4.105 + RETURNS ecircle 4.106 + LANGUAGE C IMMUTABLE STRICT 4.107 + AS '$libdir/latlon-v0008', 'pgl_ecircle_recv'; 4.108 + 4.109 +CREATE FUNCTION epoint_send(epoint) 4.110 + RETURNS bytea 4.111 + LANGUAGE C IMMUTABLE STRICT 4.112 + AS '$libdir/latlon-v0008', 'pgl_epoint_send'; 4.113 + 4.114 +CREATE FUNCTION ebox_send(ebox) 4.115 + RETURNS bytea 4.116 + LANGUAGE C IMMUTABLE STRICT 4.117 + AS '$libdir/latlon-v0008', 'pgl_ebox_send'; 4.118 + 4.119 +CREATE FUNCTION ecircle_send(ecircle) 4.120 + RETURNS bytea 4.121 + LANGUAGE C IMMUTABLE STRICT 4.122 + AS '$libdir/latlon-v0008', 'pgl_ecircle_send'; 4.123 + 4.124 + 4.125 +----------------------------------------------- 4.126 +-- type definitions of dummy index key types -- 4.127 +----------------------------------------------- 4.128 + 4.129 +CREATE TYPE ekey_point ( 4.130 + internallength = 8, 4.131 + input = ekey_point_in_dummy, 4.132 + output = ekey_point_out_dummy, 4.133 + alignment = char ); 4.134 + 4.135 +CREATE TYPE ekey_area ( 4.136 + internallength = 9, 4.137 + input = ekey_area_in_dummy, 4.138 + output = ekey_area_out_dummy, 4.139 + alignment = char ); 4.140 + 4.141 + 4.142 +------------------------------------------ 4.143 +-- definitions of geographic data types -- 4.144 +------------------------------------------ 4.145 + 4.146 +CREATE TYPE epoint ( 4.147 + internallength = 16, 4.148 + input = epoint_in, 4.149 + output = epoint_out, 4.150 + receive = epoint_recv, 4.151 + send = epoint_send, 4.152 + alignment = double ); 4.153 + 4.154 +CREATE TYPE ebox ( 4.155 + internallength = 32, 4.156 + input = ebox_in, 4.157 + output = ebox_out, 4.158 + receive = ebox_recv, 4.159 + send = ebox_send, 4.160 + alignment = double ); 4.161 + 4.162 +CREATE TYPE ecircle ( 4.163 + internallength = 24, 4.164 + input = ecircle_in, 4.165 + output = ecircle_out, 4.166 + receive = ecircle_recv, 4.167 + send = ecircle_send, 4.168 + alignment = double ); 4.169 + 4.170 +CREATE TYPE ecluster ( 4.171 + internallength = VARIABLE, 4.172 + input = ecluster_in, 4.173 + output = ecluster_out, 4.174 + alignment = double, 4.175 + storage = external ); 4.176 + 4.177 + 4.178 +-------------------- 4.179 +-- B-tree support -- 4.180 +-------------------- 4.181 + 4.182 +-- begin of B-tree support for epoint 4.183 + 4.184 +CREATE FUNCTION epoint_btree_lt(epoint, epoint) 4.185 + RETURNS boolean 4.186 + LANGUAGE C IMMUTABLE STRICT 4.187 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_lt'; 4.188 + 4.189 +CREATE FUNCTION epoint_btree_le(epoint, epoint) 4.190 + RETURNS boolean 4.191 + LANGUAGE C IMMUTABLE STRICT 4.192 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_le'; 4.193 + 4.194 +CREATE FUNCTION epoint_btree_eq(epoint, epoint) 4.195 + RETURNS boolean 4.196 + LANGUAGE C IMMUTABLE STRICT 4.197 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_eq'; 4.198 + 4.199 +CREATE FUNCTION epoint_btree_ne(epoint, epoint) 4.200 + RETURNS boolean 4.201 + LANGUAGE C IMMUTABLE STRICT 4.202 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ne'; 4.203 + 4.204 +CREATE FUNCTION epoint_btree_ge(epoint, epoint) 4.205 + RETURNS boolean 4.206 + LANGUAGE C IMMUTABLE STRICT 4.207 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ge'; 4.208 + 4.209 +CREATE FUNCTION epoint_btree_gt(epoint, epoint) 4.210 + RETURNS boolean 4.211 + LANGUAGE C IMMUTABLE STRICT 4.212 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_gt'; 4.213 + 4.214 +CREATE OPERATOR <<< ( 4.215 + leftarg = epoint, 4.216 + rightarg = epoint, 4.217 + procedure = epoint_btree_lt, 4.218 + commutator = >>>, 4.219 + negator = >>>=, 4.220 + restrict = scalarltsel, 4.221 + join = scalarltjoinsel 4.222 +); 4.223 + 4.224 +CREATE OPERATOR <<<= ( 4.225 + leftarg = epoint, 4.226 + rightarg = epoint, 4.227 + procedure = epoint_btree_le, 4.228 + commutator = >>>=, 4.229 + negator = >>>, 4.230 + restrict = scalarltsel, 4.231 + join = scalarltjoinsel 4.232 +); 4.233 + 4.234 +CREATE OPERATOR = ( 4.235 + leftarg = epoint, 4.236 + rightarg = epoint, 4.237 + procedure = epoint_btree_eq, 4.238 + commutator = =, 4.239 + negator = <>, 4.240 + restrict = eqsel, 4.241 + join = eqjoinsel, 4.242 + merges 4.243 +); 4.244 + 4.245 +CREATE OPERATOR <> ( 4.246 + leftarg = epoint, 4.247 + rightarg = epoint, 4.248 + procedure = epoint_btree_eq, 4.249 + commutator = <>, 4.250 + negator = =, 4.251 + restrict = neqsel, 4.252 + join = neqjoinsel 4.253 +); 4.254 + 4.255 +CREATE OPERATOR >>>= ( 4.256 + leftarg = epoint, 4.257 + rightarg = epoint, 4.258 + procedure = epoint_btree_ge, 4.259 + commutator = <<<=, 4.260 + negator = <<<, 4.261 + restrict = scalargtsel, 4.262 + join = scalargtjoinsel 4.263 +); 4.264 + 4.265 +CREATE OPERATOR >>> ( 4.266 + leftarg = epoint, 4.267 + rightarg = epoint, 4.268 + procedure = epoint_btree_gt, 4.269 + commutator = <<<, 4.270 + negator = <<<=, 4.271 + restrict = scalargtsel, 4.272 + join = scalargtjoinsel 4.273 +); 4.274 + 4.275 +CREATE FUNCTION epoint_btree_cmp(epoint, epoint) 4.276 + RETURNS int4 4.277 + LANGUAGE C IMMUTABLE STRICT 4.278 + AS '$libdir/latlon-v0008', 'pgl_btree_epoint_cmp'; 4.279 + 4.280 +CREATE OPERATOR CLASS epoint_btree_ops 4.281 + DEFAULT FOR TYPE epoint USING btree AS 4.282 + OPERATOR 1 <<< , 4.283 + OPERATOR 2 <<<= , 4.284 + OPERATOR 3 = , 4.285 + OPERATOR 4 >>>= , 4.286 + OPERATOR 5 >>> , 4.287 + FUNCTION 1 epoint_btree_cmp(epoint, epoint); 4.288 + 4.289 +-- end of B-tree support for epoint 4.290 + 4.291 +-- begin of B-tree support for ebox 4.292 + 4.293 +CREATE FUNCTION ebox_btree_lt(ebox, ebox) 4.294 + RETURNS boolean 4.295 + LANGUAGE C IMMUTABLE STRICT 4.296 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_lt'; 4.297 + 4.298 +CREATE FUNCTION ebox_btree_le(ebox, ebox) 4.299 + RETURNS boolean 4.300 + LANGUAGE C IMMUTABLE STRICT 4.301 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_le'; 4.302 + 4.303 +CREATE FUNCTION ebox_btree_eq(ebox, ebox) 4.304 + RETURNS boolean 4.305 + LANGUAGE C IMMUTABLE STRICT 4.306 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_eq'; 4.307 + 4.308 +CREATE FUNCTION ebox_btree_ne(ebox, ebox) 4.309 + RETURNS boolean 4.310 + LANGUAGE C IMMUTABLE STRICT 4.311 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ne'; 4.312 + 4.313 +CREATE FUNCTION ebox_btree_ge(ebox, ebox) 4.314 + RETURNS boolean 4.315 + LANGUAGE C IMMUTABLE STRICT 4.316 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ge'; 4.317 + 4.318 +CREATE FUNCTION ebox_btree_gt(ebox, ebox) 4.319 + RETURNS boolean 4.320 + LANGUAGE C IMMUTABLE STRICT 4.321 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_gt'; 4.322 + 4.323 +CREATE OPERATOR <<< ( 4.324 + leftarg = ebox, 4.325 + rightarg = ebox, 4.326 + procedure = ebox_btree_lt, 4.327 + commutator = >>>, 4.328 + negator = >>>=, 4.329 + restrict = scalarltsel, 4.330 + join = scalarltjoinsel 4.331 +); 4.332 + 4.333 +CREATE OPERATOR <<<= ( 4.334 + leftarg = ebox, 4.335 + rightarg = ebox, 4.336 + procedure = ebox_btree_le, 4.337 + commutator = >>>=, 4.338 + negator = >>>, 4.339 + restrict = scalarltsel, 4.340 + join = scalarltjoinsel 4.341 +); 4.342 + 4.343 +CREATE OPERATOR = ( 4.344 + leftarg = ebox, 4.345 + rightarg = ebox, 4.346 + procedure = ebox_btree_eq, 4.347 + commutator = =, 4.348 + negator = <>, 4.349 + restrict = eqsel, 4.350 + join = eqjoinsel, 4.351 + merges 4.352 +); 4.353 + 4.354 +CREATE OPERATOR <> ( 4.355 + leftarg = ebox, 4.356 + rightarg = ebox, 4.357 + procedure = ebox_btree_eq, 4.358 + commutator = <>, 4.359 + negator = =, 4.360 + restrict = neqsel, 4.361 + join = neqjoinsel 4.362 +); 4.363 + 4.364 +CREATE OPERATOR >>>= ( 4.365 + leftarg = ebox, 4.366 + rightarg = ebox, 4.367 + procedure = ebox_btree_ge, 4.368 + commutator = <<<=, 4.369 + negator = <<<, 4.370 + restrict = scalargtsel, 4.371 + join = scalargtjoinsel 4.372 +); 4.373 + 4.374 +CREATE OPERATOR >>> ( 4.375 + leftarg = ebox, 4.376 + rightarg = ebox, 4.377 + procedure = ebox_btree_gt, 4.378 + commutator = <<<, 4.379 + negator = <<<=, 4.380 + restrict = scalargtsel, 4.381 + join = scalargtjoinsel 4.382 +); 4.383 + 4.384 +CREATE FUNCTION ebox_btree_cmp(ebox, ebox) 4.385 + RETURNS int4 4.386 + LANGUAGE C IMMUTABLE STRICT 4.387 + AS '$libdir/latlon-v0008', 'pgl_btree_ebox_cmp'; 4.388 + 4.389 +CREATE OPERATOR CLASS ebox_btree_ops 4.390 + DEFAULT FOR TYPE ebox USING btree AS 4.391 + OPERATOR 1 <<< , 4.392 + OPERATOR 2 <<<= , 4.393 + OPERATOR 3 = , 4.394 + OPERATOR 4 >>>= , 4.395 + OPERATOR 5 >>> , 4.396 + FUNCTION 1 ebox_btree_cmp(ebox, ebox); 4.397 + 4.398 +-- end of B-tree support for ebox 4.399 + 4.400 +-- begin of B-tree support for ecircle 4.401 + 4.402 +CREATE FUNCTION ecircle_btree_lt(ecircle, ecircle) 4.403 + RETURNS boolean 4.404 + LANGUAGE C IMMUTABLE STRICT 4.405 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_lt'; 4.406 + 4.407 +CREATE FUNCTION ecircle_btree_le(ecircle, ecircle) 4.408 + RETURNS boolean 4.409 + LANGUAGE C IMMUTABLE STRICT 4.410 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_le'; 4.411 + 4.412 +CREATE FUNCTION ecircle_btree_eq(ecircle, ecircle) 4.413 + RETURNS boolean 4.414 + LANGUAGE C IMMUTABLE STRICT 4.415 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_eq'; 4.416 + 4.417 +CREATE FUNCTION ecircle_btree_ne(ecircle, ecircle) 4.418 + RETURNS boolean 4.419 + LANGUAGE C IMMUTABLE STRICT 4.420 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ne'; 4.421 + 4.422 +CREATE FUNCTION ecircle_btree_ge(ecircle, ecircle) 4.423 + RETURNS boolean 4.424 + LANGUAGE C IMMUTABLE STRICT 4.425 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ge'; 4.426 + 4.427 +CREATE FUNCTION ecircle_btree_gt(ecircle, ecircle) 4.428 + RETURNS boolean 4.429 + LANGUAGE C IMMUTABLE STRICT 4.430 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_gt'; 4.431 + 4.432 +CREATE OPERATOR <<< ( 4.433 + leftarg = ecircle, 4.434 + rightarg = ecircle, 4.435 + procedure = ecircle_btree_lt, 4.436 + commutator = >>>, 4.437 + negator = >>>=, 4.438 + restrict = scalarltsel, 4.439 + join = scalarltjoinsel 4.440 +); 4.441 + 4.442 +CREATE OPERATOR <<<= ( 4.443 + leftarg = ecircle, 4.444 + rightarg = ecircle, 4.445 + procedure = ecircle_btree_le, 4.446 + commutator = >>>=, 4.447 + negator = >>>, 4.448 + restrict = scalarltsel, 4.449 + join = scalarltjoinsel 4.450 +); 4.451 + 4.452 +CREATE OPERATOR = ( 4.453 + leftarg = ecircle, 4.454 + rightarg = ecircle, 4.455 + procedure = ecircle_btree_eq, 4.456 + commutator = =, 4.457 + negator = <>, 4.458 + restrict = eqsel, 4.459 + join = eqjoinsel, 4.460 + merges 4.461 +); 4.462 + 4.463 +CREATE OPERATOR <> ( 4.464 + leftarg = ecircle, 4.465 + rightarg = ecircle, 4.466 + procedure = ecircle_btree_eq, 4.467 + commutator = <>, 4.468 + negator = =, 4.469 + restrict = neqsel, 4.470 + join = neqjoinsel 4.471 +); 4.472 + 4.473 +CREATE OPERATOR >>>= ( 4.474 + leftarg = ecircle, 4.475 + rightarg = ecircle, 4.476 + procedure = ecircle_btree_ge, 4.477 + commutator = <<<=, 4.478 + negator = <<<, 4.479 + restrict = scalargtsel, 4.480 + join = scalargtjoinsel 4.481 +); 4.482 + 4.483 +CREATE OPERATOR >>> ( 4.484 + leftarg = ecircle, 4.485 + rightarg = ecircle, 4.486 + procedure = ecircle_btree_gt, 4.487 + commutator = <<<, 4.488 + negator = <<<=, 4.489 + restrict = scalargtsel, 4.490 + join = scalargtjoinsel 4.491 +); 4.492 + 4.493 +CREATE FUNCTION ecircle_btree_cmp(ecircle, ecircle) 4.494 + RETURNS int4 4.495 + LANGUAGE C IMMUTABLE STRICT 4.496 + AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_cmp'; 4.497 + 4.498 +CREATE OPERATOR CLASS ecircle_btree_ops 4.499 + DEFAULT FOR TYPE ecircle USING btree AS 4.500 + OPERATOR 1 <<< , 4.501 + OPERATOR 2 <<<= , 4.502 + OPERATOR 3 = , 4.503 + OPERATOR 4 >>>= , 4.504 + OPERATOR 5 >>> , 4.505 + FUNCTION 1 ecircle_btree_cmp(ecircle, ecircle); 4.506 + 4.507 +-- end of B-tree support for ecircle 4.508 + 4.509 + 4.510 +---------------- 4.511 +-- type casts -- 4.512 +---------------- 4.513 + 4.514 +CREATE FUNCTION cast_epoint_to_ebox(epoint) 4.515 + RETURNS ebox 4.516 + LANGUAGE C IMMUTABLE STRICT 4.517 + AS '$libdir/latlon-v0008', 'pgl_epoint_to_ebox'; 4.518 + 4.519 +CREATE CAST (epoint AS ebox) WITH FUNCTION cast_epoint_to_ebox(epoint); 4.520 + 4.521 +CREATE FUNCTION cast_epoint_to_ecircle(epoint) 4.522 + RETURNS ecircle 4.523 + LANGUAGE C IMMUTABLE STRICT 4.524 + AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecircle'; 4.525 + 4.526 +CREATE CAST (epoint AS ecircle) WITH FUNCTION cast_epoint_to_ecircle(epoint); 4.527 + 4.528 +CREATE FUNCTION cast_epoint_to_ecluster(epoint) 4.529 + RETURNS ecluster 4.530 + LANGUAGE C IMMUTABLE STRICT 4.531 + AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecluster'; 4.532 + 4.533 +CREATE CAST (epoint AS ecluster) WITH FUNCTION cast_epoint_to_ecluster(epoint); 4.534 + 4.535 +CREATE FUNCTION cast_ebox_to_ecluster(ebox) 4.536 + RETURNS ecluster 4.537 + LANGUAGE C IMMUTABLE STRICT 4.538 + AS '$libdir/latlon-v0008', 'pgl_ebox_to_ecluster'; 4.539 + 4.540 +CREATE CAST (ebox AS ecluster) WITH FUNCTION cast_ebox_to_ecluster(ebox); 4.541 + 4.542 + 4.543 +--------------------------- 4.544 +-- constructor functions -- 4.545 +--------------------------- 4.546 + 4.547 +CREATE FUNCTION epoint(float8, float8) 4.548 + RETURNS epoint 4.549 + LANGUAGE C IMMUTABLE STRICT 4.550 + AS '$libdir/latlon-v0008', 'pgl_create_epoint'; 4.551 + 4.552 +CREATE FUNCTION epoint_latlon(float8, float8) 4.553 + RETURNS epoint 4.554 + LANGUAGE SQL IMMUTABLE STRICT AS $$ 4.555 + SELECT epoint($1, $2) 4.556 + $$; 4.557 + 4.558 +CREATE FUNCTION epoint_lonlat(float8, float8) 4.559 + RETURNS epoint 4.560 + LANGUAGE SQL IMMUTABLE STRICT AS $$ 4.561 + SELECT epoint($2, $1) 4.562 + $$; 4.563 + 4.564 +CREATE FUNCTION empty_ebox() 4.565 + RETURNS ebox 4.566 + LANGUAGE C IMMUTABLE STRICT 4.567 + AS '$libdir/latlon-v0008', 'pgl_create_empty_ebox'; 4.568 + 4.569 +CREATE FUNCTION ebox(float8, float8, float8, float8) 4.570 + RETURNS ebox 4.571 + LANGUAGE C IMMUTABLE STRICT 4.572 + AS '$libdir/latlon-v0008', 'pgl_create_ebox'; 4.573 + 4.574 +CREATE FUNCTION ebox(epoint, epoint) 4.575 + RETURNS ebox 4.576 + LANGUAGE C IMMUTABLE STRICT 4.577 + AS '$libdir/latlon-v0008', 'pgl_create_ebox_from_epoints'; 4.578 + 4.579 +CREATE FUNCTION ecircle(float8, float8, float8) 4.580 + RETURNS ecircle 4.581 + LANGUAGE C IMMUTABLE STRICT 4.582 + AS '$libdir/latlon-v0008', 'pgl_create_ecircle'; 4.583 + 4.584 +CREATE FUNCTION ecircle(epoint, float8) 4.585 + RETURNS ecircle 4.586 + LANGUAGE C IMMUTABLE STRICT 4.587 + AS '$libdir/latlon-v0008', 'pgl_create_ecircle_from_epoint'; 4.588 + 4.589 +CREATE FUNCTION ecluster_concat(ecluster[]) 4.590 + RETURNS ecluster 4.591 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.592 + SELECT array_to_string($1, ' ')::ecluster 4.593 + $$; 4.594 + 4.595 +CREATE FUNCTION ecluster_concat(ecluster, ecluster) 4.596 + RETURNS ecluster 4.597 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.598 + SELECT ($1::text || ' ' || $2::text)::ecluster 4.599 + $$; 4.600 + 4.601 +CREATE FUNCTION ecluster_create_multipoint(epoint[]) 4.602 + RETURNS ecluster 4.603 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.604 + SELECT 4.605 + array_to_string(array_agg('point (' || unnest || ')'), ' ')::ecluster 4.606 + FROM unnest($1) 4.607 + $$; 4.608 + 4.609 +CREATE FUNCTION ecluster_create_path(epoint[]) 4.610 + RETURNS ecluster 4.611 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.612 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 4.613 + ('path (' || array_to_string($1, ' ') || ')')::ecluster 4.614 + END 4.615 + FROM array_to_string($1, ' ') AS "str" 4.616 + $$; 4.617 + 4.618 +CREATE FUNCTION ecluster_create_outline(epoint[]) 4.619 + RETURNS ecluster 4.620 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.621 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 4.622 + ('outline (' || array_to_string($1, ' ') || ')')::ecluster 4.623 + END 4.624 + FROM array_to_string($1, ' ') AS "str" 4.625 + $$; 4.626 + 4.627 +CREATE FUNCTION ecluster_create_polygon(epoint[]) 4.628 + RETURNS ecluster 4.629 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.630 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 4.631 + ('polygon (' || array_to_string($1, ' ') || ')')::ecluster 4.632 + END 4.633 + FROM array_to_string($1, ' ') AS "str" 4.634 + $$; 4.635 + 4.636 + 4.637 +---------------------- 4.638 +-- getter functions -- 4.639 +---------------------- 4.640 + 4.641 +CREATE FUNCTION latitude(epoint) 4.642 + RETURNS float8 4.643 + LANGUAGE C IMMUTABLE STRICT 4.644 + AS '$libdir/latlon-v0008', 'pgl_epoint_lat'; 4.645 + 4.646 +CREATE FUNCTION longitude(epoint) 4.647 + RETURNS float8 4.648 + LANGUAGE C IMMUTABLE STRICT 4.649 + AS '$libdir/latlon-v0008', 'pgl_epoint_lon'; 4.650 + 4.651 +CREATE FUNCTION min_latitude(ebox) 4.652 + RETURNS float8 4.653 + LANGUAGE C IMMUTABLE STRICT 4.654 + AS '$libdir/latlon-v0008', 'pgl_ebox_lat_min'; 4.655 + 4.656 +CREATE FUNCTION max_latitude(ebox) 4.657 + RETURNS float8 4.658 + LANGUAGE C IMMUTABLE STRICT 4.659 + AS '$libdir/latlon-v0008', 'pgl_ebox_lat_max'; 4.660 + 4.661 +CREATE FUNCTION min_longitude(ebox) 4.662 + RETURNS float8 4.663 + LANGUAGE C IMMUTABLE STRICT 4.664 + AS '$libdir/latlon-v0008', 'pgl_ebox_lon_min'; 4.665 + 4.666 +CREATE FUNCTION max_longitude(ebox) 4.667 + RETURNS float8 4.668 + LANGUAGE C IMMUTABLE STRICT 4.669 + AS '$libdir/latlon-v0008', 'pgl_ebox_lon_max'; 4.670 + 4.671 +CREATE FUNCTION center(ecircle) 4.672 + RETURNS epoint 4.673 + LANGUAGE C IMMUTABLE STRICT 4.674 + AS '$libdir/latlon-v0008', 'pgl_ecircle_center'; 4.675 + 4.676 +CREATE FUNCTION radius(ecircle) 4.677 + RETURNS float8 4.678 + LANGUAGE C IMMUTABLE STRICT 4.679 + AS '$libdir/latlon-v0008', 'pgl_ecircle_radius'; 4.680 + 4.681 +CREATE FUNCTION ecluster_extract_points(ecluster) 4.682 + RETURNS SETOF epoint 4.683 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.684 + SELECT "match"[2]::epoint 4.685 + FROM regexp_matches($1::text, e'(^| )point \\(([^)]+)\\)', 'g') AS "match" 4.686 + $$; 4.687 + 4.688 +CREATE FUNCTION ecluster_extract_paths(ecluster) 4.689 + RETURNS SETOF epoint[] 4.690 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.691 + SELECT ( 4.692 + SELECT array_agg("m2"[1]::epoint) 4.693 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 4.694 + ) 4.695 + FROM regexp_matches($1::text, e'(^| )path \\(([^)]+)\\)', 'g') AS "m1" 4.696 + $$; 4.697 + 4.698 +CREATE FUNCTION ecluster_extract_outlines(ecluster) 4.699 + RETURNS SETOF epoint[] 4.700 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.701 + SELECT ( 4.702 + SELECT array_agg("m2"[1]::epoint) 4.703 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 4.704 + ) 4.705 + FROM regexp_matches($1::text, e'(^| )outline \\(([^)]+)\\)', 'g') AS "m1" 4.706 + $$; 4.707 + 4.708 +CREATE FUNCTION ecluster_extract_polygons(ecluster) 4.709 + RETURNS SETOF epoint[] 4.710 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.711 + SELECT ( 4.712 + SELECT array_agg("m2"[1]::epoint) 4.713 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 4.714 + ) 4.715 + FROM regexp_matches($1::text, e'(^| )polygon \\(([^)]+)\\)', 'g') AS "m1" 4.716 + $$; 4.717 + 4.718 + 4.719 +--------------- 4.720 +-- operators -- 4.721 +--------------- 4.722 + 4.723 +CREATE FUNCTION epoint_ebox_overlap_proc(epoint, ebox) 4.724 + RETURNS boolean 4.725 + LANGUAGE C IMMUTABLE STRICT 4.726 + AS '$libdir/latlon-v0008', 'pgl_epoint_ebox_overlap'; 4.727 + 4.728 +CREATE FUNCTION epoint_ecircle_overlap_proc(epoint, ecircle) 4.729 + RETURNS boolean 4.730 + LANGUAGE C IMMUTABLE STRICT 4.731 + AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_overlap'; 4.732 + 4.733 +CREATE FUNCTION epoint_ecluster_overlap_proc(epoint, ecluster) 4.734 + RETURNS boolean 4.735 + LANGUAGE C IMMUTABLE STRICT 4.736 + AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_overlap'; 4.737 + 4.738 +CREATE FUNCTION epoint_ecluster_may_overlap_proc(epoint, ecluster) 4.739 + RETURNS boolean 4.740 + LANGUAGE C IMMUTABLE STRICT 4.741 + AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_may_overlap'; 4.742 + 4.743 +CREATE FUNCTION ebox_overlap_proc(ebox, ebox) 4.744 + RETURNS boolean 4.745 + LANGUAGE C IMMUTABLE STRICT 4.746 + AS '$libdir/latlon-v0008', 'pgl_ebox_overlap'; 4.747 + 4.748 +CREATE FUNCTION ebox_ecircle_may_overlap_proc(ebox, ecircle) 4.749 + RETURNS boolean 4.750 + LANGUAGE C IMMUTABLE STRICT 4.751 + AS '$libdir/latlon-v0008', 'pgl_ebox_ecircle_may_overlap'; 4.752 + 4.753 +CREATE FUNCTION ebox_ecluster_may_overlap_proc(ebox, ecluster) 4.754 + RETURNS boolean 4.755 + LANGUAGE C IMMUTABLE STRICT 4.756 + AS '$libdir/latlon-v0008', 'pgl_ebox_ecluster_may_overlap'; 4.757 + 4.758 +CREATE FUNCTION ecircle_overlap_proc(ecircle, ecircle) 4.759 + RETURNS boolean 4.760 + LANGUAGE C IMMUTABLE STRICT 4.761 + AS '$libdir/latlon-v0008', 'pgl_ecircle_overlap'; 4.762 + 4.763 +CREATE FUNCTION ecircle_ecluster_overlap_proc(ecircle, ecluster) 4.764 + RETURNS boolean 4.765 + LANGUAGE C IMMUTABLE STRICT 4.766 + AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_overlap'; 4.767 + 4.768 +CREATE FUNCTION ecircle_ecluster_may_overlap_proc(ecircle, ecluster) 4.769 + RETURNS boolean 4.770 + LANGUAGE C IMMUTABLE STRICT 4.771 + AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_may_overlap'; 4.772 + 4.773 +CREATE FUNCTION ecluster_overlap_proc(ecluster, ecluster) 4.774 + RETURNS boolean 4.775 + LANGUAGE C IMMUTABLE STRICT 4.776 + AS '$libdir/latlon-v0008', 'pgl_ecluster_overlap'; 4.777 + 4.778 +CREATE FUNCTION ecluster_may_overlap_proc(ecluster, ecluster) 4.779 + RETURNS boolean 4.780 + LANGUAGE C IMMUTABLE STRICT 4.781 + AS '$libdir/latlon-v0008', 'pgl_ecluster_may_overlap'; 4.782 + 4.783 +CREATE FUNCTION ecluster_contains_proc(ecluster, ecluster) 4.784 + RETURNS boolean 4.785 + LANGUAGE C IMMUTABLE STRICT 4.786 + AS '$libdir/latlon-v0008', 'pgl_ecluster_contains'; 4.787 + 4.788 +CREATE FUNCTION epoint_distance_proc(epoint, epoint) 4.789 + RETURNS float8 4.790 + LANGUAGE C IMMUTABLE STRICT 4.791 + AS '$libdir/latlon-v0008', 'pgl_epoint_distance'; 4.792 + 4.793 +CREATE FUNCTION epoint_ecircle_distance_proc(epoint, ecircle) 4.794 + RETURNS float8 4.795 + LANGUAGE C IMMUTABLE STRICT 4.796 + AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_distance'; 4.797 + 4.798 +CREATE FUNCTION epoint_ecluster_distance_proc(epoint, ecluster) 4.799 + RETURNS float8 4.800 + LANGUAGE C IMMUTABLE STRICT 4.801 + AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_distance'; 4.802 + 4.803 +CREATE FUNCTION ecircle_distance_proc(ecircle, ecircle) 4.804 + RETURNS float8 4.805 + LANGUAGE C IMMUTABLE STRICT 4.806 + AS '$libdir/latlon-v0008', 'pgl_ecircle_distance'; 4.807 + 4.808 +CREATE FUNCTION ecircle_ecluster_distance_proc(ecircle, ecluster) 4.809 + RETURNS float8 4.810 + LANGUAGE C IMMUTABLE STRICT 4.811 + AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_distance'; 4.812 + 4.813 +CREATE FUNCTION ecluster_distance_proc(ecluster, ecluster) 4.814 + RETURNS float8 4.815 + LANGUAGE C IMMUTABLE STRICT 4.816 + AS '$libdir/latlon-v0008', 'pgl_ecluster_distance'; 4.817 + 4.818 +CREATE OPERATOR && ( 4.819 + leftarg = epoint, 4.820 + rightarg = ebox, 4.821 + procedure = epoint_ebox_overlap_proc, 4.822 + commutator = &&, 4.823 + restrict = areasel, 4.824 + join = areajoinsel 4.825 +); 4.826 + 4.827 +CREATE FUNCTION epoint_ebox_overlap_commutator(ebox, epoint) 4.828 + RETURNS boolean 4.829 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 4.830 + 4.831 +CREATE OPERATOR && ( 4.832 + leftarg = ebox, 4.833 + rightarg = epoint, 4.834 + procedure = epoint_ebox_overlap_commutator, 4.835 + commutator = &&, 4.836 + restrict = areasel, 4.837 + join = areajoinsel 4.838 +); 4.839 + 4.840 +CREATE OPERATOR && ( 4.841 + leftarg = epoint, 4.842 + rightarg = ecircle, 4.843 + procedure = epoint_ecircle_overlap_proc, 4.844 + commutator = &&, 4.845 + restrict = areasel, 4.846 + join = areajoinsel 4.847 +); 4.848 + 4.849 +CREATE FUNCTION epoint_ecircle_overlap_commutator(ecircle, epoint) 4.850 + RETURNS boolean 4.851 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 4.852 + 4.853 +CREATE OPERATOR && ( 4.854 + leftarg = ecircle, 4.855 + rightarg = epoint, 4.856 + procedure = epoint_ecircle_overlap_commutator, 4.857 + commutator = &&, 4.858 + restrict = areasel, 4.859 + join = areajoinsel 4.860 +); 4.861 + 4.862 +CREATE OPERATOR && ( 4.863 + leftarg = epoint, 4.864 + rightarg = ecluster, 4.865 + procedure = epoint_ecluster_overlap_proc, 4.866 + commutator = &&, 4.867 + restrict = areasel, 4.868 + join = areajoinsel 4.869 +); 4.870 + 4.871 +CREATE FUNCTION epoint_ecluster_overlap_commutator(ecluster, epoint) 4.872 + RETURNS boolean 4.873 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 4.874 + 4.875 +CREATE OPERATOR && ( 4.876 + leftarg = ecluster, 4.877 + rightarg = epoint, 4.878 + procedure = epoint_ecluster_overlap_commutator, 4.879 + commutator = &&, 4.880 + restrict = areasel, 4.881 + join = areajoinsel 4.882 +); 4.883 + 4.884 +CREATE OPERATOR && ( 4.885 + leftarg = ebox, 4.886 + rightarg = ebox, 4.887 + procedure = ebox_overlap_proc, 4.888 + commutator = &&, 4.889 + restrict = areasel, 4.890 + join = areajoinsel 4.891 +); 4.892 + 4.893 +CREATE OPERATOR && ( 4.894 + leftarg = ecircle, 4.895 + rightarg = ecircle, 4.896 + procedure = ecircle_overlap_proc, 4.897 + commutator = &&, 4.898 + restrict = areasel, 4.899 + join = areajoinsel 4.900 +); 4.901 + 4.902 +CREATE OPERATOR && ( 4.903 + leftarg = ecircle, 4.904 + rightarg = ecluster, 4.905 + procedure = ecircle_ecluster_overlap_proc, 4.906 + commutator = &&, 4.907 + restrict = areasel, 4.908 + join = areajoinsel 4.909 +); 4.910 + 4.911 +CREATE FUNCTION ecircle_ecluster_overlap_commutator(ecluster, ecircle) 4.912 + RETURNS boolean 4.913 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 4.914 + 4.915 +CREATE OPERATOR && ( 4.916 + leftarg = ecluster, 4.917 + rightarg = ecircle, 4.918 + procedure = ecircle_ecluster_overlap_commutator, 4.919 + commutator = &&, 4.920 + restrict = areasel, 4.921 + join = areajoinsel 4.922 +); 4.923 + 4.924 +CREATE OPERATOR && ( 4.925 + leftarg = ecluster, 4.926 + rightarg = ecluster, 4.927 + procedure = ecluster_overlap_proc, 4.928 + commutator = &&, 4.929 + restrict = areasel, 4.930 + join = areajoinsel 4.931 +); 4.932 + 4.933 +CREATE FUNCTION ebox_ecircle_overlap_castwrap(ebox, ecircle) 4.934 + RETURNS boolean 4.935 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2'; 4.936 + 4.937 +CREATE OPERATOR && ( 4.938 + leftarg = ebox, 4.939 + rightarg = ecircle, 4.940 + procedure = ebox_ecircle_overlap_castwrap, 4.941 + commutator = &&, 4.942 + restrict = areasel, 4.943 + join = areajoinsel 4.944 +); 4.945 + 4.946 +CREATE FUNCTION ebox_ecircle_overlap_castwrap(ecircle, ebox) 4.947 + RETURNS boolean 4.948 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster'; 4.949 + 4.950 +CREATE OPERATOR && ( 4.951 + leftarg = ecircle, 4.952 + rightarg = ebox, 4.953 + procedure = ebox_ecircle_overlap_castwrap, 4.954 + commutator = &&, 4.955 + restrict = areasel, 4.956 + join = areajoinsel 4.957 +); 4.958 + 4.959 +CREATE FUNCTION ebox_ecluster_overlap_castwrap(ebox, ecluster) 4.960 + RETURNS boolean 4.961 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2'; 4.962 + 4.963 +CREATE OPERATOR && ( 4.964 + leftarg = ebox, 4.965 + rightarg = ecluster, 4.966 + procedure = ebox_ecluster_overlap_castwrap, 4.967 + commutator = &&, 4.968 + restrict = areasel, 4.969 + join = areajoinsel 4.970 +); 4.971 + 4.972 +CREATE FUNCTION ebox_ecluster_overlap_castwrap(ecluster, ebox) 4.973 + RETURNS boolean 4.974 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster'; 4.975 + 4.976 +CREATE OPERATOR && ( 4.977 + leftarg = ecluster, 4.978 + rightarg = ebox, 4.979 + procedure = ebox_ecluster_overlap_castwrap, 4.980 + commutator = &&, 4.981 + restrict = areasel, 4.982 + join = areajoinsel 4.983 +); 4.984 + 4.985 +CREATE OPERATOR &&+ ( 4.986 + leftarg = epoint, 4.987 + rightarg = ecluster, 4.988 + procedure = epoint_ecluster_may_overlap_proc, 4.989 + commutator = &&+, 4.990 + restrict = areasel, 4.991 + join = areajoinsel 4.992 +); 4.993 + 4.994 +CREATE FUNCTION epoint_ecluster_may_overlap_commutator(ecluster, epoint) 4.995 + RETURNS boolean 4.996 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1'; 4.997 + 4.998 +CREATE OPERATOR &&+ ( 4.999 + leftarg = ecluster, 4.1000 + rightarg = epoint, 4.1001 + procedure = epoint_ecluster_may_overlap_commutator, 4.1002 + commutator = &&+, 4.1003 + restrict = areasel, 4.1004 + join = areajoinsel 4.1005 +); 4.1006 + 4.1007 +CREATE OPERATOR &&+ ( 4.1008 + leftarg = ebox, 4.1009 + rightarg = ecircle, 4.1010 + procedure = ebox_ecircle_may_overlap_proc, 4.1011 + commutator = &&+, 4.1012 + restrict = areasel, 4.1013 + join = areajoinsel 4.1014 +); 4.1015 + 4.1016 +CREATE FUNCTION ebox_ecircle_may_overlap_commutator(ecircle, ebox) 4.1017 + RETURNS boolean 4.1018 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1'; 4.1019 + 4.1020 +CREATE OPERATOR &&+ ( 4.1021 + leftarg = ecircle, 4.1022 + rightarg = ebox, 4.1023 + procedure = ebox_ecircle_may_overlap_commutator, 4.1024 + commutator = &&+, 4.1025 + restrict = areasel, 4.1026 + join = areajoinsel 4.1027 +); 4.1028 + 4.1029 +CREATE OPERATOR &&+ ( 4.1030 + leftarg = ebox, 4.1031 + rightarg = ecluster, 4.1032 + procedure = ebox_ecluster_may_overlap_proc, 4.1033 + commutator = &&+, 4.1034 + restrict = areasel, 4.1035 + join = areajoinsel 4.1036 +); 4.1037 + 4.1038 +CREATE FUNCTION ebox_ecluster_may_overlap_commutator(ecluster, ebox) 4.1039 + RETURNS boolean 4.1040 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1'; 4.1041 + 4.1042 +CREATE OPERATOR &&+ ( 4.1043 + leftarg = ecluster, 4.1044 + rightarg = ebox, 4.1045 + procedure = ebox_ecluster_may_overlap_commutator, 4.1046 + commutator = &&+, 4.1047 + restrict = areasel, 4.1048 + join = areajoinsel 4.1049 +); 4.1050 + 4.1051 +CREATE OPERATOR &&+ ( 4.1052 + leftarg = ecircle, 4.1053 + rightarg = ecluster, 4.1054 + procedure = ecircle_ecluster_may_overlap_proc, 4.1055 + commutator = &&+, 4.1056 + restrict = areasel, 4.1057 + join = areajoinsel 4.1058 +); 4.1059 + 4.1060 +CREATE FUNCTION ecircle_ecluster_may_overlap_commutator(ecluster, ecircle) 4.1061 + RETURNS boolean 4.1062 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1'; 4.1063 + 4.1064 +CREATE OPERATOR &&+ ( 4.1065 + leftarg = ecluster, 4.1066 + rightarg = ecircle, 4.1067 + procedure = ecircle_ecluster_may_overlap_commutator, 4.1068 + commutator = &&+, 4.1069 + restrict = areasel, 4.1070 + join = areajoinsel 4.1071 +); 4.1072 + 4.1073 +CREATE OPERATOR &&+ ( 4.1074 + leftarg = ecluster, 4.1075 + rightarg = ecluster, 4.1076 + procedure = ecluster_may_overlap_proc, 4.1077 + commutator = &&+, 4.1078 + restrict = areasel, 4.1079 + join = areajoinsel 4.1080 +); 4.1081 + 4.1082 +CREATE OPERATOR @> ( 4.1083 + leftarg = ebox, 4.1084 + rightarg = epoint, 4.1085 + procedure = epoint_ebox_overlap_commutator, 4.1086 + commutator = <@, 4.1087 + restrict = areasel, 4.1088 + join = areajoinsel 4.1089 +); 4.1090 + 4.1091 +CREATE OPERATOR <@ ( 4.1092 + leftarg = epoint, 4.1093 + rightarg = ebox, 4.1094 + procedure = epoint_ebox_overlap_proc, 4.1095 + commutator = @>, 4.1096 + restrict = areasel, 4.1097 + join = areajoinsel 4.1098 +); 4.1099 + 4.1100 +CREATE OPERATOR @> ( 4.1101 + leftarg = ecluster, 4.1102 + rightarg = epoint, 4.1103 + procedure = epoint_ecluster_overlap_commutator, 4.1104 + commutator = <@, 4.1105 + restrict = areasel, 4.1106 + join = areajoinsel 4.1107 +); 4.1108 + 4.1109 +CREATE OPERATOR <@ ( 4.1110 + leftarg = epoint, 4.1111 + rightarg = ecluster, 4.1112 + procedure = epoint_ecluster_overlap_proc, 4.1113 + commutator = <@, 4.1114 + restrict = areasel, 4.1115 + join = areajoinsel 4.1116 +); 4.1117 + 4.1118 +CREATE OPERATOR @> ( 4.1119 + leftarg = ecluster, 4.1120 + rightarg = ecluster, 4.1121 + procedure = ecluster_contains_proc, 4.1122 + commutator = <@, 4.1123 + restrict = areasel, 4.1124 + join = areajoinsel 4.1125 +); 4.1126 + 4.1127 +CREATE FUNCTION ecluster_contains_commutator(ecluster, ecluster) 4.1128 + RETURNS boolean 4.1129 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1'; 4.1130 + 4.1131 +CREATE OPERATOR <@ ( 4.1132 + leftarg = ecluster, 4.1133 + rightarg = ecluster, 4.1134 + procedure = ecluster_contains_commutator, 4.1135 + commutator = @>, 4.1136 + restrict = areasel, 4.1137 + join = areajoinsel 4.1138 +); 4.1139 + 4.1140 +CREATE FUNCTION ebox_contains_castwrap(ebox, ebox) 4.1141 + RETURNS boolean 4.1142 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2::ecluster'; 4.1143 + 4.1144 +CREATE OPERATOR @> ( 4.1145 + leftarg = ebox, 4.1146 + rightarg = ebox, 4.1147 + procedure = ebox_contains_castwrap, 4.1148 + commutator = <@, 4.1149 + restrict = areasel, 4.1150 + join = areajoinsel 4.1151 +); 4.1152 + 4.1153 +CREATE FUNCTION ebox_contains_swapped_castwrap(ebox, ebox) 4.1154 + RETURNS boolean 4.1155 + LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1::ecluster'; 4.1156 + 4.1157 +CREATE OPERATOR <@ ( 4.1158 + leftarg = ebox, 4.1159 + rightarg = ebox, 4.1160 + procedure = ebox_contains_swapped_castwrap, 4.1161 + commutator = @>, 4.1162 + restrict = areasel, 4.1163 + join = areajoinsel 4.1164 +); 4.1165 + 4.1166 +CREATE FUNCTION ebox_ecluster_contains_castwrap(ebox, ecluster) 4.1167 + RETURNS boolean 4.1168 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2'; 4.1169 + 4.1170 +CREATE OPERATOR @> ( 4.1171 + leftarg = ebox, 4.1172 + rightarg = ecluster, 4.1173 + procedure = ebox_ecluster_contains_castwrap, 4.1174 + commutator = <@, 4.1175 + restrict = areasel, 4.1176 + join = areajoinsel 4.1177 +); 4.1178 + 4.1179 +CREATE FUNCTION ebox_ecluster_contains_castwrap(ecluster, ebox) 4.1180 + RETURNS boolean 4.1181 + LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1'; 4.1182 + 4.1183 +CREATE OPERATOR <@ ( 4.1184 + leftarg = ecluster, 4.1185 + rightarg = ebox, 4.1186 + procedure = ebox_ecluster_contains_castwrap, 4.1187 + commutator = @>, 4.1188 + restrict = areasel, 4.1189 + join = areajoinsel 4.1190 +); 4.1191 + 4.1192 +CREATE FUNCTION ecluster_ebox_contains_castwrap(ecluster, ebox) 4.1193 + RETURNS boolean 4.1194 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 @> $2::ecluster'; 4.1195 + 4.1196 +CREATE OPERATOR @> ( 4.1197 + leftarg = ecluster, 4.1198 + rightarg = ebox, 4.1199 + procedure = ecluster_ebox_contains_castwrap, 4.1200 + commutator = <@, 4.1201 + restrict = areasel, 4.1202 + join = areajoinsel 4.1203 +); 4.1204 + 4.1205 +CREATE FUNCTION ecluster_ebox_contains_castwrap(ebox, ecluster) 4.1206 + RETURNS boolean 4.1207 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1::ecluster'; 4.1208 + 4.1209 +CREATE OPERATOR <@ ( 4.1210 + leftarg = ebox, 4.1211 + rightarg = ecluster, 4.1212 + procedure = ecluster_ebox_contains_castwrap, 4.1213 + commutator = @>, 4.1214 + restrict = areasel, 4.1215 + join = areajoinsel 4.1216 +); 4.1217 + 4.1218 +CREATE OPERATOR <-> ( 4.1219 + leftarg = epoint, 4.1220 + rightarg = epoint, 4.1221 + procedure = epoint_distance_proc, 4.1222 + commutator = <-> 4.1223 +); 4.1224 + 4.1225 +CREATE OPERATOR <-> ( 4.1226 + leftarg = epoint, 4.1227 + rightarg = ecircle, 4.1228 + procedure = epoint_ecircle_distance_proc, 4.1229 + commutator = <-> 4.1230 +); 4.1231 + 4.1232 +CREATE FUNCTION epoint_ecircle_distance_commutator(ecircle, epoint) 4.1233 + RETURNS float8 4.1234 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 4.1235 + 4.1236 +CREATE OPERATOR <-> ( 4.1237 + leftarg = ecircle, 4.1238 + rightarg = epoint, 4.1239 + procedure = epoint_ecircle_distance_commutator, 4.1240 + commutator = <-> 4.1241 +); 4.1242 + 4.1243 +CREATE OPERATOR <-> ( 4.1244 + leftarg = epoint, 4.1245 + rightarg = ecluster, 4.1246 + procedure = epoint_ecluster_distance_proc, 4.1247 + commutator = <-> 4.1248 +); 4.1249 + 4.1250 +CREATE FUNCTION epoint_ecluster_distance_commutator(ecluster, epoint) 4.1251 + RETURNS float8 4.1252 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 4.1253 + 4.1254 +CREATE OPERATOR <-> ( 4.1255 + leftarg = ecluster, 4.1256 + rightarg = epoint, 4.1257 + procedure = epoint_ecluster_distance_commutator, 4.1258 + commutator = <-> 4.1259 +); 4.1260 + 4.1261 +CREATE OPERATOR <-> ( 4.1262 + leftarg = ecircle, 4.1263 + rightarg = ecircle, 4.1264 + procedure = ecircle_distance_proc, 4.1265 + commutator = <-> 4.1266 +); 4.1267 + 4.1268 +CREATE OPERATOR <-> ( 4.1269 + leftarg = ecircle, 4.1270 + rightarg = ecluster, 4.1271 + procedure = ecircle_ecluster_distance_proc, 4.1272 + commutator = <-> 4.1273 +); 4.1274 + 4.1275 +CREATE FUNCTION ecircle_ecluster_distance_commutator(ecluster, ecircle) 4.1276 + RETURNS float8 4.1277 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 4.1278 + 4.1279 +CREATE OPERATOR <-> ( 4.1280 + leftarg = ecluster, 4.1281 + rightarg = ecircle, 4.1282 + procedure = ecircle_ecluster_distance_commutator, 4.1283 + commutator = <-> 4.1284 +); 4.1285 + 4.1286 +CREATE OPERATOR <-> ( 4.1287 + leftarg = ecluster, 4.1288 + rightarg = ecluster, 4.1289 + procedure = ecluster_distance_proc, 4.1290 + commutator = <-> 4.1291 +); 4.1292 + 4.1293 +CREATE FUNCTION epoint_ebox_distance_castwrap(epoint, ebox) 4.1294 + RETURNS float8 4.1295 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster'; 4.1296 + 4.1297 +CREATE OPERATOR <-> ( 4.1298 + leftarg = epoint, 4.1299 + rightarg = ebox, 4.1300 + procedure = epoint_ebox_distance_castwrap, 4.1301 + commutator = <-> 4.1302 +); 4.1303 + 4.1304 +CREATE FUNCTION epoint_ebox_distance_castwrap(ebox, epoint) 4.1305 + RETURNS float8 4.1306 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2'; 4.1307 + 4.1308 +CREATE OPERATOR <-> ( 4.1309 + leftarg = ebox, 4.1310 + rightarg = epoint, 4.1311 + procedure = epoint_ebox_distance_castwrap, 4.1312 + commutator = <-> 4.1313 +); 4.1314 + 4.1315 +CREATE FUNCTION ebox_distance_castwrap(ebox, ebox) 4.1316 + RETURNS float8 4.1317 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2::ecluster'; 4.1318 + 4.1319 +CREATE OPERATOR <-> ( 4.1320 + leftarg = ebox, 4.1321 + rightarg = ebox, 4.1322 + procedure = ebox_distance_castwrap, 4.1323 + commutator = <-> 4.1324 +); 4.1325 + 4.1326 +CREATE FUNCTION ebox_ecircle_distance_castwrap(ebox, ecircle) 4.1327 + RETURNS float8 4.1328 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2'; 4.1329 + 4.1330 +CREATE OPERATOR <-> ( 4.1331 + leftarg = ebox, 4.1332 + rightarg = ecircle, 4.1333 + procedure = ebox_ecircle_distance_castwrap, 4.1334 + commutator = <-> 4.1335 +); 4.1336 + 4.1337 +CREATE FUNCTION ebox_ecircle_distance_castwrap(ecircle, ebox) 4.1338 + RETURNS float8 4.1339 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster'; 4.1340 + 4.1341 +CREATE OPERATOR <-> ( 4.1342 + leftarg = ecircle, 4.1343 + rightarg = ebox, 4.1344 + procedure = ebox_ecircle_distance_castwrap, 4.1345 + commutator = <-> 4.1346 +); 4.1347 + 4.1348 +CREATE FUNCTION ebox_ecluster_distance_castwrap(ebox, ecluster) 4.1349 + RETURNS float8 4.1350 + LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2'; 4.1351 + 4.1352 +CREATE OPERATOR <-> ( 4.1353 + leftarg = ebox, 4.1354 + rightarg = ecluster, 4.1355 + procedure = ebox_ecluster_distance_castwrap, 4.1356 + commutator = <-> 4.1357 +); 4.1358 + 4.1359 +CREATE FUNCTION ebox_ecluster_distance_castwrap(ecluster, ebox) 4.1360 + RETURNS float8 4.1361 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster'; 4.1362 + 4.1363 +CREATE OPERATOR <-> ( 4.1364 + leftarg = ecluster, 4.1365 + rightarg = ebox, 4.1366 + procedure = ebox_ecluster_distance_castwrap, 4.1367 + commutator = <-> 4.1368 +); 4.1369 + 4.1370 + 4.1371 +--------------------- 4.1372 +-- other functions -- 4.1373 +--------------------- 4.1374 + 4.1375 +CREATE FUNCTION monte_carlo_area(ecluster, int4 = 10000) 4.1376 + RETURNS float8 4.1377 + LANGUAGE C STRICT 4.1378 + AS '$libdir/latlon-v0008', 'pgl_ecluster_monte_carlo_area'; 4.1379 + 4.1380 +CREATE FUNCTION fair_distance(ecluster, epoint, int4 = 10000) 4.1381 + RETURNS float8 4.1382 + LANGUAGE C STRICT 4.1383 + AS '$libdir/latlon-v0008', 'pgl_ecluster_epoint_fair_distance'; 4.1384 + 4.1385 + 4.1386 +---------------- 4.1387 +-- GiST index -- 4.1388 +---------------- 4.1389 + 4.1390 +CREATE FUNCTION pgl_gist_consistent(internal, internal, smallint, oid, internal) 4.1391 + RETURNS boolean 4.1392 + LANGUAGE C STRICT 4.1393 + AS '$libdir/latlon-v0008', 'pgl_gist_consistent'; 4.1394 + 4.1395 +CREATE FUNCTION pgl_gist_union(internal, internal) 4.1396 + RETURNS internal 4.1397 + LANGUAGE C STRICT 4.1398 + AS '$libdir/latlon-v0008', 'pgl_gist_union'; 4.1399 + 4.1400 +CREATE FUNCTION pgl_gist_compress_epoint(internal) 4.1401 + RETURNS internal 4.1402 + LANGUAGE C STRICT 4.1403 + AS '$libdir/latlon-v0008', 'pgl_gist_compress_epoint'; 4.1404 + 4.1405 +CREATE FUNCTION pgl_gist_compress_ecircle(internal) 4.1406 + RETURNS internal 4.1407 + LANGUAGE C STRICT 4.1408 + AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecircle'; 4.1409 + 4.1410 +CREATE FUNCTION pgl_gist_compress_ecluster(internal) 4.1411 + RETURNS internal 4.1412 + LANGUAGE C STRICT 4.1413 + AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecluster'; 4.1414 + 4.1415 +CREATE FUNCTION pgl_gist_decompress(internal) 4.1416 + RETURNS internal 4.1417 + LANGUAGE C STRICT 4.1418 + AS '$libdir/latlon-v0008', 'pgl_gist_decompress'; 4.1419 + 4.1420 +CREATE FUNCTION pgl_gist_penalty(internal, internal, internal) 4.1421 + RETURNS internal 4.1422 + LANGUAGE C STRICT 4.1423 + AS '$libdir/latlon-v0008', 'pgl_gist_penalty'; 4.1424 + 4.1425 +CREATE FUNCTION pgl_gist_picksplit(internal, internal) 4.1426 + RETURNS internal 4.1427 + LANGUAGE C STRICT 4.1428 + AS '$libdir/latlon-v0008', 'pgl_gist_picksplit'; 4.1429 + 4.1430 +CREATE FUNCTION pgl_gist_same(internal, internal, internal) 4.1431 + RETURNS internal 4.1432 + LANGUAGE C STRICT 4.1433 + AS '$libdir/latlon-v0008', 'pgl_gist_same'; 4.1434 + 4.1435 +CREATE FUNCTION pgl_gist_distance(internal, internal, smallint, oid) 4.1436 + RETURNS internal 4.1437 + LANGUAGE C STRICT 4.1438 + AS '$libdir/latlon-v0008', 'pgl_gist_distance'; 4.1439 + 4.1440 +CREATE OPERATOR CLASS epoint_ops 4.1441 + DEFAULT FOR TYPE epoint USING gist AS 4.1442 + OPERATOR 11 = , 4.1443 + OPERATOR 22 && (epoint, ebox), 4.1444 + OPERATOR 222 <@ (epoint, ebox), 4.1445 + OPERATOR 23 && (epoint, ecircle), 4.1446 + OPERATOR 24 && (epoint, ecluster), 4.1447 + OPERATOR 124 &&+ (epoint, ecluster), 4.1448 + OPERATOR 224 <@ (epoint, ecluster), 4.1449 + OPERATOR 31 <-> (epoint, epoint) FOR ORDER BY float_ops, 4.1450 + OPERATOR 32 <-> (epoint, ebox) FOR ORDER BY float_ops, 4.1451 + OPERATOR 33 <-> (epoint, ecircle) FOR ORDER BY float_ops, 4.1452 + OPERATOR 34 <-> (epoint, ecluster) FOR ORDER BY float_ops, 4.1453 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 4.1454 + FUNCTION 2 pgl_gist_union(internal, internal), 4.1455 + FUNCTION 3 pgl_gist_compress_epoint(internal), 4.1456 + FUNCTION 4 pgl_gist_decompress(internal), 4.1457 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 4.1458 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 4.1459 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 4.1460 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 4.1461 + STORAGE ekey_point; 4.1462 + 4.1463 +CREATE OPERATOR CLASS ecircle_ops 4.1464 + DEFAULT FOR TYPE ecircle USING gist AS 4.1465 + OPERATOR 13 = , 4.1466 + OPERATOR 21 && (ecircle, epoint), 4.1467 + OPERATOR 22 && (ecircle, ebox), 4.1468 + OPERATOR 122 &&+ (ecircle, ebox), 4.1469 + OPERATOR 23 && (ecircle, ecircle), 4.1470 + OPERATOR 24 && (ecircle, ecluster), 4.1471 + OPERATOR 124 &&+ (ecircle, ecluster), 4.1472 + OPERATOR 31 <-> (ecircle, epoint) FOR ORDER BY float_ops, 4.1473 + OPERATOR 32 <-> (ecircle, ebox) FOR ORDER BY float_ops, 4.1474 + OPERATOR 33 <-> (ecircle, ecircle) FOR ORDER BY float_ops, 4.1475 + OPERATOR 34 <-> (ecircle, ecluster) FOR ORDER BY float_ops, 4.1476 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 4.1477 + FUNCTION 2 pgl_gist_union(internal, internal), 4.1478 + FUNCTION 3 pgl_gist_compress_ecircle(internal), 4.1479 + FUNCTION 4 pgl_gist_decompress(internal), 4.1480 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 4.1481 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 4.1482 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 4.1483 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 4.1484 + STORAGE ekey_area; 4.1485 + 4.1486 +CREATE OPERATOR CLASS ecluster_ops 4.1487 + DEFAULT FOR TYPE ecluster USING gist AS 4.1488 + OPERATOR 21 && (ecluster, epoint), 4.1489 + OPERATOR 121 &&+ (ecluster, epoint), 4.1490 + OPERATOR 221 @> (ecluster, epoint), 4.1491 + OPERATOR 22 && (ecluster, ebox), 4.1492 + OPERATOR 122 &&+ (ecluster, ebox), 4.1493 + OPERATOR 222 @> (ecluster, ebox), 4.1494 + OPERATOR 322 <@ (ecluster, ebox), 4.1495 + OPERATOR 23 && (ecluster, ecircle), 4.1496 + OPERATOR 123 &&+ (ecluster, ecircle), 4.1497 + OPERATOR 24 && (ecluster, ecluster), 4.1498 + OPERATOR 124 &&+ (ecluster, ecluster), 4.1499 + OPERATOR 224 @> (ecluster, ecluster), 4.1500 + OPERATOR 324 <@ (ecluster, ecluster), 4.1501 + OPERATOR 31 <-> (ecluster, epoint) FOR ORDER BY float_ops, 4.1502 + OPERATOR 32 <-> (ecluster, ebox) FOR ORDER BY float_ops, 4.1503 + OPERATOR 33 <-> (ecluster, ecircle) FOR ORDER BY float_ops, 4.1504 + OPERATOR 34 <-> (ecluster, ecluster) FOR ORDER BY float_ops, 4.1505 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 4.1506 + FUNCTION 2 pgl_gist_union(internal, internal), 4.1507 + FUNCTION 3 pgl_gist_compress_ecluster(internal), 4.1508 + FUNCTION 4 pgl_gist_decompress(internal), 4.1509 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 4.1510 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 4.1511 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 4.1512 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 4.1513 + STORAGE ekey_area; 4.1514 + 4.1515 + 4.1516 +--------------------- 4.1517 +-- alias functions -- 4.1518 +--------------------- 4.1519 + 4.1520 +CREATE FUNCTION distance(epoint, epoint) 4.1521 + RETURNS float8 4.1522 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2'; 4.1523 + 4.1524 +CREATE FUNCTION distance(ecluster, epoint) 4.1525 + RETURNS float8 4.1526 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2'; 4.1527 + 4.1528 +CREATE FUNCTION distance_within(epoint, epoint, float8) 4.1529 + RETURNS boolean 4.1530 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)'; 4.1531 + 4.1532 +CREATE FUNCTION distance_within(ecluster, epoint, float8) 4.1533 + RETURNS boolean 4.1534 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)'; 4.1535 + 4.1536 + 4.1537 +-------------------------------- 4.1538 +-- other data storage formats -- 4.1539 +-------------------------------- 4.1540 + 4.1541 +CREATE FUNCTION coords_to_epoint(float8, float8, text = 'epoint') 4.1542 + RETURNS epoint 4.1543 + LANGUAGE plpgsql IMMUTABLE STRICT AS $$ 4.1544 + DECLARE 4.1545 + "result" epoint; 4.1546 + BEGIN 4.1547 + IF $3 = 'epoint_lonlat' THEN 4.1548 + -- avoid dynamic command execution for better performance 4.1549 + RETURN epoint($2, $1); 4.1550 + END IF; 4.1551 + IF $3 = 'epoint' OR $3 = 'epoint_latlon' THEN 4.1552 + -- avoid dynamic command execution for better performance 4.1553 + RETURN epoint($1, $2); 4.1554 + END IF; 4.1555 + EXECUTE 'SELECT ' || $3 || '($1, $2)' INTO STRICT "result" USING $1, $2; 4.1556 + RETURN "result"; 4.1557 + END; 4.1558 + $$; 4.1559 + 4.1560 +CREATE FUNCTION GeoJSON_LinearRing_vertices(jsonb, text = 'epoint_lonlat') 4.1561 + RETURNS SETOF jsonb 4.1562 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.1563 + SELECT "result" FROM 4.1564 + ( SELECT jsonb_array_length($1) - 1 ) AS "lastindex_row" ("lastindex") 4.1565 + CROSS JOIN LATERAL jsonb_array_elements( 4.1566 + CASE WHEN 4.1567 + coords_to_epoint( 4.1568 + ($1->0->>0)::float8, 4.1569 + ($1->0->>1)::float8, 4.1570 + $2 4.1571 + ) = coords_to_epoint( 4.1572 + ($1->"lastindex"->>0)::float8, 4.1573 + ($1->"lastindex"->>1)::float8, 4.1574 + $2 4.1575 + ) 4.1576 + THEN 4.1577 + $1 - "lastindex" 4.1578 + ELSE 4.1579 + $1 4.1580 + END 4.1581 + ) AS "result_row" ("result") 4.1582 + $$; 4.1583 + 4.1584 +CREATE FUNCTION GeoJSON_to_epoint(jsonb, text = 'epoint_lonlat') 4.1585 + RETURNS epoint 4.1586 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.1587 + SELECT CASE 4.1588 + WHEN $1->>'type' = 'Point' THEN 4.1589 + coords_to_epoint( 4.1590 + ($1->'coordinates'->>0)::float8, 4.1591 + ($1->'coordinates'->>1)::float8, 4.1592 + $2 4.1593 + ) 4.1594 + WHEN $1->>'type' = 'Feature' THEN 4.1595 + GeoJSON_to_epoint($1->'geometry', $2) 4.1596 + ELSE 4.1597 + NULL 4.1598 + END 4.1599 + $$; 4.1600 + 4.1601 +CREATE FUNCTION GeoJSON_to_ecluster(jsonb, text = 'epoint_lonlat') 4.1602 + RETURNS ecluster 4.1603 + LANGUAGE sql IMMUTABLE STRICT AS $$ 4.1604 + SELECT CASE $1->>'type' 4.1605 + WHEN 'Point' THEN 4.1606 + coords_to_epoint( 4.1607 + ($1->'coordinates'->>0)::float8, 4.1608 + ($1->'coordinates'->>1)::float8, 4.1609 + $2 4.1610 + )::ecluster 4.1611 + WHEN 'MultiPoint' THEN 4.1612 + ( SELECT ecluster_create_multipoint(array_agg( 4.1613 + coords_to_epoint( 4.1614 + ("coord"->>0)::float8, 4.1615 + ("coord"->>1)::float8, 4.1616 + $2 4.1617 + ) 4.1618 + )) 4.1619 + FROM jsonb_array_elements($1->'coordinates') AS "coord" 4.1620 + ) 4.1621 + WHEN 'LineString' THEN 4.1622 + ( SELECT ecluster_create_path(array_agg( 4.1623 + coords_to_epoint( 4.1624 + ("coord"->>0)::float8, 4.1625 + ("coord"->>1)::float8, 4.1626 + $2 4.1627 + ) 4.1628 + )) 4.1629 + FROM jsonb_array_elements($1->'coordinates') AS "coord" 4.1630 + ) 4.1631 + WHEN 'MultiLineString' THEN 4.1632 + ( SELECT ecluster_concat(array_agg( 4.1633 + ( SELECT ecluster_create_path(array_agg( 4.1634 + coords_to_epoint( 4.1635 + ("coord"->>0)::float8, 4.1636 + ("coord"->>1)::float8, 4.1637 + $2 4.1638 + ) 4.1639 + )) 4.1640 + FROM jsonb_array_elements("coord_array") AS "coord" 4.1641 + ) 4.1642 + )) 4.1643 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array" 4.1644 + ) 4.1645 + WHEN 'Polygon' THEN 4.1646 + ( SELECT ecluster_concat(array_agg( 4.1647 + ( SELECT ecluster_create_polygon(array_agg( 4.1648 + coords_to_epoint( 4.1649 + ("coord"->>0)::float8, 4.1650 + ("coord"->>1)::float8, 4.1651 + $2 4.1652 + ) 4.1653 + )) 4.1654 + FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord" 4.1655 + ) 4.1656 + )) 4.1657 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array" 4.1658 + ) 4.1659 + WHEN 'MultiPolygon' THEN 4.1660 + ( SELECT ecluster_concat(array_agg( 4.1661 + ( SELECT ecluster_concat(array_agg( 4.1662 + ( SELECT ecluster_create_polygon(array_agg( 4.1663 + coords_to_epoint( 4.1664 + ("coord"->>0)::float8, 4.1665 + ("coord"->>1)::float8, 4.1666 + $2 4.1667 + ) 4.1668 + )) 4.1669 + FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord" 4.1670 + ) 4.1671 + )) 4.1672 + FROM jsonb_array_elements("coord_array_array") AS "coord_array" 4.1673 + ) 4.1674 + )) 4.1675 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array_array" 4.1676 + ) 4.1677 + WHEN 'GeometryCollection' THEN 4.1678 + ( SELECT ecluster_concat(array_agg( 4.1679 + GeoJSON_to_ecluster("geometry", $2) 4.1680 + )) 4.1681 + FROM jsonb_array_elements($1->'geometries') AS "geometry" 4.1682 + ) 4.1683 + WHEN 'Feature' THEN 4.1684 + GeoJSON_to_ecluster($1->'geometry', $2) 4.1685 + WHEN 'FeatureCollection' THEN 4.1686 + ( SELECT ecluster_concat(array_agg( 4.1687 + GeoJSON_to_ecluster("feature", $2) 4.1688 + )) 4.1689 + FROM jsonb_array_elements($1->'features') AS "feature" 4.1690 + ) 4.1691 + ELSE 4.1692 + NULL 4.1693 + END 4.1694 + $$; 4.1695 +
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 5.2 +++ b/latlon-v0008.c Tue Oct 25 18:44:43 2016 +0200 5.3 @@ -0,0 +1,3285 @@ 5.4 + 5.5 +/*-------------* 5.6 + * C prelude * 5.7 + *-------------*/ 5.8 + 5.9 +#include "postgres.h" 5.10 +#include "fmgr.h" 5.11 +#include "libpq/pqformat.h" 5.12 +#include "access/gist.h" 5.13 +#include "access/stratnum.h" 5.14 +#include "utils/array.h" 5.15 +#include <math.h> 5.16 + 5.17 +#ifdef PG_MODULE_MAGIC 5.18 +PG_MODULE_MAGIC; 5.19 +#endif 5.20 + 5.21 +#if INT_MAX < 2147483647 5.22 +#error Expected int type to be at least 32 bit wide 5.23 +#endif 5.24 + 5.25 + 5.26 +/*---------------------------------* 5.27 + * distance calculation on earth * 5.28 + * (using WGS-84 spheroid) * 5.29 + *---------------------------------*/ 5.30 + 5.31 +/* WGS-84 spheroid with following parameters: 5.32 + semi-major axis a = 6378137 5.33 + semi-minor axis b = a * (1 - 1/298.257223563) 5.34 + estimated diameter = 2 * (2*a+b)/3 5.35 +*/ 5.36 +#define PGL_SPHEROID_A 6378137.0 /* semi major axis */ 5.37 +#define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */ 5.38 +#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F)) 5.39 +#define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \ 5.40 + PGL_SPHEROID_B * PGL_SPHEROID_B ) / \ 5.41 + ( PGL_SPHEROID_A * PGL_SPHEROID_A ) ) 5.42 +#define PGL_SUBEPS2 (1.0-PGL_EPS2) 5.43 +#define PGL_RADIUS ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0) 5.44 +#define PGL_DIAMETER (2.0 * PGL_RADIUS) 5.45 +#define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */ 5.46 +#define PGL_MAXDIST (PGL_RADIUS * M_PI) /* maximum distance */ 5.47 +#define PGL_FADELIMIT (PGL_MAXDIST / 3.0) /* 1/6 circumference */ 5.48 + 5.49 +/* calculate distance between two points on earth (given in degrees) */ 5.50 +static inline double pgl_distance( 5.51 + double lat1, double lon1, double lat2, double lon2 5.52 +) { 5.53 + float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin; 5.54 + float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t; 5.55 + /* normalize delta longitude (lon2 > 0 && lon1 = 0) */ 5.56 + /* lon1 = 0 (not used anymore) */ 5.57 + lon2 = fabs(lon2-lon1); 5.58 + /* convert to radians (first divide, then multiply) */ 5.59 + lat1 = (lat1 / 180.0) * M_PI; 5.60 + lat2 = (lat2 / 180.0) * M_PI; 5.61 + lon2 = (lon2 / 180.0) * M_PI; 5.62 + /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point 5.63 + operations (lon2 >= lon1 is already ensured in a previous step) */ 5.64 + if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; } 5.65 + /* calculate 3d coordinates on scaled ellipsoid which has an average diameter 5.66 + of 1.0 */ 5.67 + lat1cos = cos(lat1); lat1sin = sin(lat1); 5.68 + lat2cos = cos(lat2); lat2sin = sin(lat2); 5.69 + lon2cos = cos(lon2); lon2sin = sin(lon2); 5.70 + nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin); 5.71 + nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin); 5.72 + x1 = nphi1 * lat1cos; 5.73 + z1 = nphi1 * PGL_SUBEPS2 * lat1sin; 5.74 + x2 = nphi2 * lat2cos * lon2cos; 5.75 + y2 = nphi2 * lat2cos * lon2sin; 5.76 + z2 = nphi2 * PGL_SUBEPS2 * lat2sin; 5.77 + /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */ 5.78 + g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1)); 5.79 + /* convert tunnel distance through scaled ellipsoid to approximated surface 5.80 + distance on original ellipsoid */ 5.81 + if (g > 1.0) g = 1.0; 5.82 + s = PGL_DIAMETER * asin(g); 5.83 + /* return result only if small enough to be precise (less than 1/3 of 5.84 + maximum possible distance) */ 5.85 + if (s <= PGL_FADELIMIT) return s; 5.86 + /* calculate tunnel distance to antipodal point through scaled ellipsoid */ 5.87 + g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1)); 5.88 + /* convert tunnel distance to antipodal point through scaled ellipsoid to 5.89 + approximated surface distance to antipodal point on original ellipsoid */ 5.90 + if (g > 1.0) g = 1.0; 5.91 + t = PGL_DIAMETER * asin(g); 5.92 + /* surface distance between original points can now be approximated by 5.93 + substracting antipodal distance from maximum possible distance; 5.94 + return result only if small enough (less than 1/3 of maximum possible 5.95 + distance) */ 5.96 + if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t; 5.97 + /* otherwise crossfade direct and antipodal result to ensure monotonicity */ 5.98 + return ( 5.99 + (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) / 5.100 + (s + t - 2*PGL_FADELIMIT) 5.101 + ); 5.102 +} 5.103 + 5.104 +/* finite distance that can not be reached on earth */ 5.105 +#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST) 5.106 + 5.107 + 5.108 +/*--------------------------------* 5.109 + * simple geographic data types * 5.110 + *--------------------------------*/ 5.111 + 5.112 +/* point on earth given by latitude and longitude in degrees */ 5.113 +/* (type "epoint" in SQL) */ 5.114 +typedef struct { 5.115 + double lat; /* between -90 and 90 (both inclusive) */ 5.116 + double lon; /* between -180 and 180 (both inclusive) */ 5.117 +} pgl_point; 5.118 + 5.119 +/* box delimited by two parallels and two meridians (all in degrees) */ 5.120 +/* (type "ebox" in SQL) */ 5.121 +typedef struct { 5.122 + double lat_min; /* between -90 and 90 (both inclusive) */ 5.123 + double lat_max; /* between -90 and 90 (both inclusive) */ 5.124 + double lon_min; /* between -180 and 180 (both inclusive) */ 5.125 + double lon_max; /* between -180 and 180 (both inclusive) */ 5.126 + /* if lat_min > lat_max, then box is empty */ 5.127 + /* if lon_min > lon_max, then 180th meridian is crossed */ 5.128 +} pgl_box; 5.129 + 5.130 +/* circle on earth surface (for radial searches with fixed radius) */ 5.131 +/* (type "ecircle" in SQL) */ 5.132 +typedef struct { 5.133 + pgl_point center; 5.134 + double radius; /* positive (including +0 but excluding -0), or -INFINITY */ 5.135 + /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point), 5.136 + zero radius (0) denotes a single point, 5.137 + a finite radius (0 < radius < INFINITY) denotes a filled circle, and 5.138 + a radius of INFINITY is valid and means complete coverage of earth. */ 5.139 +} pgl_circle; 5.140 + 5.141 + 5.142 +/*----------------------------------* 5.143 + * geographic "cluster" data type * 5.144 + *----------------------------------*/ 5.145 + 5.146 +/* A cluster is a collection of points, paths, outlines, and polygons. If two 5.147 + polygons in a cluster overlap, the area covered by both polygons does not 5.148 + belong to the cluster. This way, a cluster can be used to describe complex 5.149 + shapes like polygons with holes. Outlines are non-filled polygons. Paths are 5.150 + open by default (i.e. the last point in the list is not connected with the 5.151 + first point in the list). Note that each outline or polygon in a cluster 5.152 + must cover a longitude range of less than 180 degrees to avoid ambiguities. 5.153 + Areas which are larger may be split into multiple polygons. */ 5.154 + 5.155 +/* maximum number of points in a cluster */ 5.156 +/* (limited to avoid integer overflows, e.g. when allocating memory) */ 5.157 +#define PGL_CLUSTER_MAXPOINTS 16777216 5.158 + 5.159 +/* types of cluster entries */ 5.160 +#define PGL_ENTRY_POINT 1 /* a point */ 5.161 +#define PGL_ENTRY_PATH 2 /* a path from first point to last point */ 5.162 +#define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */ 5.163 +#define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */ 5.164 + 5.165 +/* Entries of a cluster are described by two different structs: pgl_newentry 5.166 + and pgl_entry. The first is used only during construction of a cluster, the 5.167 + second is used in all other cases (e.g. when reading clusters from the 5.168 + database, performing operations, etc). */ 5.169 + 5.170 +/* entry for new geographic cluster during construction of that cluster */ 5.171 +typedef struct { 5.172 + int32_t entrytype; 5.173 + int32_t npoints; 5.174 + pgl_point *points; /* pointer to an array of points (pgl_point) */ 5.175 +} pgl_newentry; 5.176 + 5.177 +/* entry of geographic cluster */ 5.178 +typedef struct { 5.179 + int32_t entrytype; /* type of entry: point, path, outline, polygon */ 5.180 + int32_t npoints; /* number of stored points (set to 1 for point entry) */ 5.181 + int32_t offset; /* offset of pgl_point array from cluster base address */ 5.182 + /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */ 5.183 +} pgl_entry; 5.184 + 5.185 +/* geographic cluster which is a collection of points, (open) paths, polygons, 5.186 + and outlines (non-filled polygons) */ 5.187 +typedef struct { 5.188 + char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */ 5.189 + int32_t nentries; /* number of stored points */ 5.190 + pgl_circle bounding; /* bounding circle */ 5.191 + /* Note: bounding circle ensures alignment of pgl_cluster for points */ 5.192 + pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */ 5.193 +} pgl_cluster; 5.194 + 5.195 +/* macro to determine memory alignment of points */ 5.196 +/* (needed to store pgl_point array after entries in pgl_cluster) */ 5.197 +typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment; 5.198 +#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned) 5.199 + 5.200 +/* macro to extract a pointer to the array of points of a cluster entry */ 5.201 +#define PGL_ENTRY_POINTS(cluster, idx) \ 5.202 + ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset)) 5.203 + 5.204 +/* convert pgl_newentry array to pgl_cluster */ 5.205 +/* NOTE: requires pgl_finalize_cluster to be called to finalize result */ 5.206 +static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) { 5.207 + int i; /* index of current entry */ 5.208 + int npoints = 0; /* number of points in whole cluster */ 5.209 + int entry_npoints; /* number of points in current entry */ 5.210 + int points_offset = PGL_POINT_ALIGNMENT * ( 5.211 + ( offsetof(pgl_cluster, entries) + 5.212 + nentries * sizeof(pgl_entry) + 5.213 + PGL_POINT_ALIGNMENT - 1 5.214 + ) / PGL_POINT_ALIGNMENT 5.215 + ); /* offset of pgl_point array from base address (considering alignment) */ 5.216 + pgl_cluster *cluster; /* new cluster to be returned */ 5.217 + /* determine total number of points */ 5.218 + for (i=0; i<nentries; i++) npoints += entries[i].npoints; 5.219 + /* allocate memory for cluster (including entries and points) */ 5.220 + cluster = palloc(points_offset + npoints * sizeof(pgl_point)); 5.221 + /* re-count total number of points to determine offset for each entry */ 5.222 + npoints = 0; 5.223 + /* copy entries and points */ 5.224 + for (i=0; i<nentries; i++) { 5.225 + /* determine number of points in entry */ 5.226 + entry_npoints = entries[i].npoints; 5.227 + /* copy entry */ 5.228 + cluster->entries[i].entrytype = entries[i].entrytype; 5.229 + cluster->entries[i].npoints = entry_npoints; 5.230 + /* calculate offset (in bytes) of pgl_point array */ 5.231 + cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point); 5.232 + /* copy points */ 5.233 + memcpy( 5.234 + PGL_ENTRY_POINTS(cluster, i), 5.235 + entries[i].points, 5.236 + entry_npoints * sizeof(pgl_point) 5.237 + ); 5.238 + /* update total number of points processed */ 5.239 + npoints += entry_npoints; 5.240 + } 5.241 + /* set number of entries in cluster */ 5.242 + cluster->nentries = nentries; 5.243 + /* set PostgreSQL header for variable sized data */ 5.244 + SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point)); 5.245 + /* return newly created cluster */ 5.246 + return cluster; 5.247 +} 5.248 + 5.249 + 5.250 +/*----------------------------------------* 5.251 + * C functions on geographic data types * 5.252 + *----------------------------------------*/ 5.253 + 5.254 +/* round latitude or longitude to 12 digits after decimal point */ 5.255 +static inline double pgl_round(double val) { 5.256 + return round(val * 1e12) / 1e12; 5.257 +} 5.258 + 5.259 +/* compare two points */ 5.260 +/* (equality when same point on earth is described, otherwise an arbitrary 5.261 + linear order) */ 5.262 +static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) { 5.263 + double lon1, lon2; /* modified longitudes for special cases */ 5.264 + /* use latitude as first ordering criterion */ 5.265 + if (point1->lat < point2->lat) return -1; 5.266 + if (point1->lat > point2->lat) return 1; 5.267 + /* determine modified longitudes (considering special case of poles and 5.268 + 180th meridian which can be described as W180 or E180) */ 5.269 + if (point1->lat == -90 || point1->lat == 90) lon1 = 0; 5.270 + else if (point1->lon == 180) lon1 = -180; 5.271 + else lon1 = point1->lon; 5.272 + if (point2->lat == -90 || point2->lat == 90) lon2 = 0; 5.273 + else if (point2->lon == 180) lon2 = -180; 5.274 + else lon2 = point2->lon; 5.275 + /* use (modified) longitude as secondary ordering criterion */ 5.276 + if (lon1 < lon2) return -1; 5.277 + if (lon1 > lon2) return 1; 5.278 + /* no difference found, points are equal */ 5.279 + return 0; 5.280 +} 5.281 + 5.282 +/* compare two boxes */ 5.283 +/* (equality when same box on earth is described, otherwise an arbitrary linear 5.284 + order) */ 5.285 +static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) { 5.286 + /* two empty boxes are equal, and an empty box is always considered "less 5.287 + than" a non-empty box */ 5.288 + if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1; 5.289 + if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0; 5.290 + if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1; 5.291 + /* use southern border as first ordering criterion */ 5.292 + if (box1->lat_min < box2->lat_min) return -1; 5.293 + if (box1->lat_min > box2->lat_min) return 1; 5.294 + /* use northern border as second ordering criterion */ 5.295 + if (box1->lat_max < box2->lat_max) return -1; 5.296 + if (box1->lat_max > box2->lat_max) return 1; 5.297 + /* use western border as third ordering criterion */ 5.298 + if (box1->lon_min < box2->lon_min) return -1; 5.299 + if (box1->lon_min > box2->lon_min) return 1; 5.300 + /* use eastern border as fourth ordering criterion */ 5.301 + if (box1->lon_max < box2->lon_max) return -1; 5.302 + if (box1->lon_max > box2->lon_max) return 1; 5.303 + /* no difference found, boxes are equal */ 5.304 + return 0; 5.305 +} 5.306 + 5.307 +/* compare two circles */ 5.308 +/* (equality when same circle on earth is described, otherwise an arbitrary 5.309 + linear order) */ 5.310 +static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) { 5.311 + /* two circles with same infinite radius (positive or negative infinity) are 5.312 + considered equal independently of center point */ 5.313 + if ( 5.314 + !isfinite(circle1->radius) && !isfinite(circle2->radius) && 5.315 + circle1->radius == circle2->radius 5.316 + ) return 0; 5.317 + /* use radius as first ordering criterion */ 5.318 + if (circle1->radius < circle2->radius) return -1; 5.319 + if (circle1->radius > circle2->radius) return 1; 5.320 + /* use center point as secondary ordering criterion */ 5.321 + return pgl_point_cmp(&(circle1->center), &(circle2->center)); 5.322 +} 5.323 + 5.324 +/* set box to empty box*/ 5.325 +static void pgl_box_set_empty(pgl_box *box) { 5.326 + box->lat_min = INFINITY; 5.327 + box->lat_max = -INFINITY; 5.328 + box->lon_min = 0; 5.329 + box->lon_max = 0; 5.330 +} 5.331 + 5.332 +/* check if point is inside a box */ 5.333 +static bool pgl_point_in_box(pgl_point *point, pgl_box *box) { 5.334 + return ( 5.335 + point->lat >= box->lat_min && point->lat <= box->lat_max && ( 5.336 + (box->lon_min > box->lon_max) ? ( 5.337 + /* box crosses 180th meridian */ 5.338 + point->lon >= box->lon_min || point->lon <= box->lon_max 5.339 + ) : ( 5.340 + /* box does not cross the 180th meridian */ 5.341 + point->lon >= box->lon_min && point->lon <= box->lon_max 5.342 + ) 5.343 + ) 5.344 + ); 5.345 +} 5.346 + 5.347 +/* check if two boxes overlap */ 5.348 +static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) { 5.349 + return ( 5.350 + box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */ 5.351 + ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) && 5.352 + ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && ( 5.353 + ( 5.354 + /* check if one and only one box crosses the 180th meridian */ 5.355 + ((box1->lon_min > box1->lon_max) ? 1 : 0) ^ 5.356 + ((box2->lon_min > box2->lon_max) ? 1 : 0) 5.357 + ) ? ( 5.358 + /* exactly one box crosses the 180th meridian */ 5.359 + box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min || 5.360 + box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max 5.361 + ) : ( 5.362 + /* no box or both boxes cross the 180th meridian */ 5.363 + ( 5.364 + (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) && 5.365 + (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max) 5.366 + ) || 5.367 + /* handle W180 == E180 */ 5.368 + ( box1->lon_min == -180 && box2->lon_max == 180 ) || 5.369 + ( box2->lon_min == -180 && box1->lon_max == 180 ) 5.370 + ) 5.371 + ) 5.372 + ); 5.373 +} 5.374 + 5.375 +/* check unambiguousness of east/west orientation of cluster entries and set 5.376 + bounding circle of cluster */ 5.377 +static bool pgl_finalize_cluster(pgl_cluster *cluster) { 5.378 + int i, j; /* i: index of entry, j: index of point in entry */ 5.379 + int npoints; /* number of points in entry */ 5.380 + int total_npoints = 0; /* total number of points in cluster */ 5.381 + pgl_point *points; /* points in entry */ 5.382 + int lon_dir; /* first point of entry west (-1) or east (+1) */ 5.383 + double lon_break = 0; /* antipodal longitude of first point in entry */ 5.384 + double lon_min, lon_max; /* covered longitude range of entry */ 5.385 + double value; /* temporary variable */ 5.386 + /* reset bounding circle center to empty circle at 0/0 coordinates */ 5.387 + cluster->bounding.center.lat = 0; 5.388 + cluster->bounding.center.lon = 0; 5.389 + cluster->bounding.radius = -INFINITY; 5.390 + /* if cluster is not empty */ 5.391 + if (cluster->nentries != 0) { 5.392 + /* iterate over all cluster entries and ensure they each cover a longitude 5.393 + range less than 180 degrees */ 5.394 + for (i=0; i<cluster->nentries; i++) { 5.395 + /* get properties of entry */ 5.396 + npoints = cluster->entries[i].npoints; 5.397 + points = PGL_ENTRY_POINTS(cluster, i); 5.398 + /* get longitude of first point of entry */ 5.399 + value = points[0].lon; 5.400 + /* initialize lon_min and lon_max with longitude of first point */ 5.401 + lon_min = value; 5.402 + lon_max = value; 5.403 + /* determine east/west orientation of first point and calculate antipodal 5.404 + longitude (Note: rounding required here) */ 5.405 + if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); } 5.406 + else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); } 5.407 + else lon_dir = 0; 5.408 + /* iterate over all other points in entry */ 5.409 + for (j=1; j<npoints; j++) { 5.410 + /* consider longitude wrap-around */ 5.411 + value = points[j].lon; 5.412 + if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360); 5.413 + else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360); 5.414 + /* update lon_min and lon_max */ 5.415 + if (value < lon_min) lon_min = value; 5.416 + else if (value > lon_max) lon_max = value; 5.417 + /* return false if 180 degrees or more are covered */ 5.418 + if (lon_max - lon_min >= 180) return false; 5.419 + } 5.420 + } 5.421 + /* iterate over all points of all entries and calculate arbitrary center 5.422 + point for bounding circle (best if center point minimizes the radius, 5.423 + but some error is allowed here) */ 5.424 + for (i=0; i<cluster->nentries; i++) { 5.425 + /* get properties of entry */ 5.426 + npoints = cluster->entries[i].npoints; 5.427 + points = PGL_ENTRY_POINTS(cluster, i); 5.428 + /* check if first entry */ 5.429 + if (i==0) { 5.430 + /* get longitude of first point of first entry in whole cluster */ 5.431 + value = points[0].lon; 5.432 + /* initialize lon_min and lon_max with longitude of first point of 5.433 + first entry in whole cluster (used to determine if whole cluster 5.434 + covers a longitude range of 180 degrees or more) */ 5.435 + lon_min = value; 5.436 + lon_max = value; 5.437 + /* determine east/west orientation of first point and calculate 5.438 + antipodal longitude (Note: rounding not necessary here) */ 5.439 + if (value < 0) { lon_dir = -1; lon_break = value + 180; } 5.440 + else if (value > 0) { lon_dir = 1; lon_break = value - 180; } 5.441 + else lon_dir = 0; 5.442 + } 5.443 + /* iterate over all points in entry */ 5.444 + for (j=0; j<npoints; j++) { 5.445 + /* longitude wrap-around (Note: rounding not necessary here) */ 5.446 + value = points[j].lon; 5.447 + if (lon_dir < 0 && value > lon_break) value -= 360; 5.448 + else if (lon_dir > 0 && value < lon_break) value += 360; 5.449 + if (value < lon_min) lon_min = value; 5.450 + else if (value > lon_max) lon_max = value; 5.451 + /* set bounding circle to cover whole earth if more than 180 degrees 5.452 + are covered */ 5.453 + if (lon_max - lon_min >= 180) { 5.454 + cluster->bounding.center.lat = 0; 5.455 + cluster->bounding.center.lon = 0; 5.456 + cluster->bounding.radius = INFINITY; 5.457 + return true; 5.458 + } 5.459 + /* add point to bounding circle center (for average calculation) */ 5.460 + cluster->bounding.center.lat += points[j].lat; 5.461 + cluster->bounding.center.lon += value; 5.462 + } 5.463 + /* count total number of points */ 5.464 + total_npoints += npoints; 5.465 + } 5.466 + /* determine average latitude and longitude of cluster */ 5.467 + cluster->bounding.center.lat /= total_npoints; 5.468 + cluster->bounding.center.lon /= total_npoints; 5.469 + /* normalize longitude of center of cluster bounding circle */ 5.470 + if (cluster->bounding.center.lon < -180) { 5.471 + cluster->bounding.center.lon += 360; 5.472 + } 5.473 + else if (cluster->bounding.center.lon > 180) { 5.474 + cluster->bounding.center.lon -= 360; 5.475 + } 5.476 + /* round bounding circle center (useful if it is used by other functions) */ 5.477 + cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat); 5.478 + cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon); 5.479 + /* calculate radius of bounding circle */ 5.480 + for (i=0; i<cluster->nentries; i++) { 5.481 + npoints = cluster->entries[i].npoints; 5.482 + points = PGL_ENTRY_POINTS(cluster, i); 5.483 + for (j=0; j<npoints; j++) { 5.484 + value = pgl_distance( 5.485 + cluster->bounding.center.lat, cluster->bounding.center.lon, 5.486 + points[j].lat, points[j].lon 5.487 + ); 5.488 + if (value > cluster->bounding.radius) cluster->bounding.radius = value; 5.489 + } 5.490 + } 5.491 + } 5.492 + /* return true (east/west orientation is unambiguous) */ 5.493 + return true; 5.494 +} 5.495 + 5.496 +/* check if point is inside cluster */ 5.497 +/* (if point is on perimeter, then true is returned if and only if 5.498 + strict == false) */ 5.499 +static bool pgl_point_in_cluster( 5.500 + pgl_point *point, 5.501 + pgl_cluster *cluster, 5.502 + bool strict 5.503 +) { 5.504 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 5.505 + int entrytype; /* type of entry */ 5.506 + int npoints; /* number of points in entry */ 5.507 + pgl_point *points; /* array of points in entry */ 5.508 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 5.509 + double lon_break = 0; /* antipodal longitude of first vertex */ 5.510 + double lat0 = point->lat; /* latitude of point */ 5.511 + double lon0; /* (adjusted) longitude of point */ 5.512 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 5.513 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 5.514 + double lon; /* longitude of intersection */ 5.515 + int counter = 0; /* counter for intersections east of point */ 5.516 + /* iterate over all entries */ 5.517 + for (i=0; i<cluster->nentries; i++) { 5.518 + /* get type of entry */ 5.519 + entrytype = cluster->entries[i].entrytype; 5.520 + /* skip all entries but polygons if perimeters are excluded */ 5.521 + if (strict && entrytype != PGL_ENTRY_POLYGON) continue; 5.522 + /* get points of entry */ 5.523 + npoints = cluster->entries[i].npoints; 5.524 + points = PGL_ENTRY_POINTS(cluster, i); 5.525 + /* determine east/west orientation of first point of entry and calculate 5.526 + antipodal longitude */ 5.527 + lon_break = points[0].lon; 5.528 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 5.529 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 5.530 + else lon_dir = 0; 5.531 + /* get longitude of point */ 5.532 + lon0 = point->lon; 5.533 + /* consider longitude wrap-around for point */ 5.534 + if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360); 5.535 + else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360); 5.536 + /* iterate over all edges and vertices */ 5.537 + for (j=0; j<npoints; j++) { 5.538 + /* return if point is on vertex of polygon */ 5.539 + if (pgl_point_cmp(point, &(points[j])) == 0) return !strict; 5.540 + /* calculate index of next vertex */ 5.541 + k = (j+1) % npoints; 5.542 + /* skip last edge unless entry is (closed) outline or polygon */ 5.543 + if ( 5.544 + k == 0 && 5.545 + entrytype != PGL_ENTRY_OUTLINE && 5.546 + entrytype != PGL_ENTRY_POLYGON 5.547 + ) continue; 5.548 + /* use previously calculated values for lat1 and lon1 if possible */ 5.549 + if (j) { 5.550 + lat1 = lat2; 5.551 + lon1 = lon2; 5.552 + } else { 5.553 + /* otherwise get latitude and longitude values of first vertex */ 5.554 + lat1 = points[0].lat; 5.555 + lon1 = points[0].lon; 5.556 + /* and consider longitude wrap-around for first vertex */ 5.557 + if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360); 5.558 + else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360); 5.559 + } 5.560 + /* get latitude and longitude of next vertex */ 5.561 + lat2 = points[k].lat; 5.562 + lon2 = points[k].lon; 5.563 + /* consider longitude wrap-around for next vertex */ 5.564 + if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360); 5.565 + else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360); 5.566 + /* return if point is on horizontal (west to east) edge of polygon */ 5.567 + if ( 5.568 + lat0 == lat1 && lat0 == lat2 && 5.569 + ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) ) 5.570 + ) return !strict; 5.571 + /* check if edge crosses east/west line of point */ 5.572 + if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) { 5.573 + /* calculate longitude of intersection */ 5.574 + lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1); 5.575 + /* return if intersection goes (approximately) through point */ 5.576 + if (pgl_round(lon) == lon0) return !strict; 5.577 + /* count intersection if east of point and entry is polygon*/ 5.578 + if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++; 5.579 + } 5.580 + } 5.581 + } 5.582 + /* return true if number of intersections is odd */ 5.583 + return counter & 1; 5.584 +} 5.585 + 5.586 +/* check if all points of the second cluster are strictly inside the first 5.587 + cluster */ 5.588 +static inline bool pgl_all_cluster_points_strictly_in_cluster( 5.589 + pgl_cluster *outer, pgl_cluster *inner 5.590 +) { 5.591 + int i, j; /* i: entry, j: point in entry */ 5.592 + int npoints; /* number of points in entry */ 5.593 + pgl_point *points; /* array of points in entry */ 5.594 + /* iterate over all entries of "inner" cluster */ 5.595 + for (i=0; i<inner->nentries; i++) { 5.596 + /* get properties of entry */ 5.597 + npoints = inner->entries[i].npoints; 5.598 + points = PGL_ENTRY_POINTS(inner, i); 5.599 + /* iterate over all points in entry of "inner" cluster */ 5.600 + for (j=0; j<npoints; j++) { 5.601 + /* return false if one point of inner cluster is not in outer cluster */ 5.602 + if (!pgl_point_in_cluster(points+j, outer, true)) return false; 5.603 + } 5.604 + } 5.605 + /* otherwise return true */ 5.606 + return true; 5.607 +} 5.608 + 5.609 +/* check if any point the second cluster is inside the first cluster */ 5.610 +static inline bool pgl_any_cluster_points_in_cluster( 5.611 + pgl_cluster *outer, pgl_cluster *inner 5.612 +) { 5.613 + int i, j; /* i: entry, j: point in entry */ 5.614 + int npoints; /* number of points in entry */ 5.615 + pgl_point *points; /* array of points in entry */ 5.616 + /* iterate over all entries of "inner" cluster */ 5.617 + for (i=0; i<inner->nentries; i++) { 5.618 + /* get properties of entry */ 5.619 + npoints = inner->entries[i].npoints; 5.620 + points = PGL_ENTRY_POINTS(inner, i); 5.621 + /* iterate over all points in entry of "inner" cluster */ 5.622 + for (j=0; j<npoints; j++) { 5.623 + /* return true if one point of inner cluster is in outer cluster */ 5.624 + if (pgl_point_in_cluster(points+j, outer, false)) return true; 5.625 + } 5.626 + } 5.627 + /* otherwise return false */ 5.628 + return false; 5.629 +} 5.630 + 5.631 +/* check if line segment strictly crosses line (not just touching) */ 5.632 +static inline bool pgl_lseg_crosses_line( 5.633 + double seg_x1, double seg_y1, double seg_x2, double seg_y2, 5.634 + double line_x1, double line_y1, double line_x2, double line_y2 5.635 +) { 5.636 + return ( 5.637 + ( 5.638 + (seg_x1-line_x1) * (line_y2-line_y1) - 5.639 + (seg_y1-line_y1) * (line_x2-line_x1) 5.640 + ) * ( 5.641 + (seg_x2-line_x1) * (line_y2-line_y1) - 5.642 + (seg_y2-line_y1) * (line_x2-line_x1) 5.643 + ) 5.644 + ) < 0; 5.645 +} 5.646 + 5.647 +/* check if paths and outlines of two clusters strictly overlap (not just 5.648 + touching) */ 5.649 +static bool pgl_outlines_overlap( 5.650 + pgl_cluster *cluster1, pgl_cluster *cluster2 5.651 +) { 5.652 + int i1, j1, k1; /* i: entry, j: point in entry, k: next point in entry */ 5.653 + int i2, j2, k2; 5.654 + int entrytype1, entrytype2; /* type of entry */ 5.655 + int npoints1, npoints2; /* number of points in entry */ 5.656 + pgl_point *points1; /* array of points in entry of cluster1 */ 5.657 + pgl_point *points2; /* array of points in entry of cluster2 */ 5.658 + int lon_dir1, lon_dir2; /* first vertex west (-1) or east (+1) */ 5.659 + double lon_break1, lon_break2; /* antipodal longitude of first vertex */ 5.660 + double lat11, lon11; /* latitude and (adjusted) longitude of vertex */ 5.661 + double lat12, lon12; /* latitude and (adjusted) longitude of next vertex */ 5.662 + double lat21, lon21; /* latitude and (adjusted) longitudes for cluster2 */ 5.663 + double lat22, lon22; 5.664 + double wrapvalue; /* temporary helper value to adjust wrap-around */ 5.665 + /* iterate over all entries of cluster1 */ 5.666 + for (i1=0; i1<cluster1->nentries; i1++) { 5.667 + /* get properties of entry in cluster1 and skip points */ 5.668 + npoints1 = cluster1->entries[i1].npoints; 5.669 + if (npoints1 < 2) continue; 5.670 + entrytype1 = cluster1->entries[i1].entrytype; 5.671 + points1 = PGL_ENTRY_POINTS(cluster1, i1); 5.672 + /* determine east/west orientation of first point and calculate antipodal 5.673 + longitude */ 5.674 + lon_break1 = points1[0].lon; 5.675 + if (lon_break1 < 0) { 5.676 + lon_dir1 = -1; 5.677 + lon_break1 = pgl_round(lon_break1 + 180); 5.678 + } else if (lon_break1 > 0) { 5.679 + lon_dir1 = 1; 5.680 + lon_break1 = pgl_round(lon_break1 - 180); 5.681 + } else lon_dir1 = 0; 5.682 + /* iterate over all edges and vertices in cluster1 */ 5.683 + for (j1=0; j1<npoints1; j1++) { 5.684 + /* calculate index of next vertex */ 5.685 + k1 = (j1+1) % npoints1; 5.686 + /* skip last edge unless entry is (closed) outline or polygon */ 5.687 + if ( 5.688 + k1 == 0 && 5.689 + entrytype1 != PGL_ENTRY_OUTLINE && 5.690 + entrytype1 != PGL_ENTRY_POLYGON 5.691 + ) continue; 5.692 + /* use previously calculated values for lat1 and lon1 if possible */ 5.693 + if (j1) { 5.694 + lat11 = lat12; 5.695 + lon11 = lon12; 5.696 + } else { 5.697 + /* otherwise get latitude and longitude values of first vertex */ 5.698 + lat11 = points1[0].lat; 5.699 + lon11 = points1[0].lon; 5.700 + /* and consider longitude wrap-around for first vertex */ 5.701 + if (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360); 5.702 + else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360); 5.703 + } 5.704 + /* get latitude and longitude of next vertex */ 5.705 + lat12 = points1[k1].lat; 5.706 + lon12 = points1[k1].lon; 5.707 + /* consider longitude wrap-around for next vertex */ 5.708 + if (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360); 5.709 + else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360); 5.710 + /* skip degenerated edges */ 5.711 + if (lat11 == lat12 && lon11 == lon12) continue; 5.712 + /* iterate over all entries of cluster2 */ 5.713 + for (i2=0; i2<cluster2->nentries; i2++) { 5.714 + /* get points and number of points of entry in cluster2 */ 5.715 + npoints2 = cluster2->entries[i2].npoints; 5.716 + if (npoints2 < 2) continue; 5.717 + entrytype2 = cluster2->entries[i2].entrytype; 5.718 + points2 = PGL_ENTRY_POINTS(cluster2, i2); 5.719 + /* determine east/west orientation of first point and calculate antipodal 5.720 + longitude */ 5.721 + lon_break2 = points2[0].lon; 5.722 + if (lon_break2 < 0) { 5.723 + lon_dir2 = -1; 5.724 + lon_break2 = pgl_round(lon_break2 + 180); 5.725 + } else if (lon_break2 > 0) { 5.726 + lon_dir2 = 1; 5.727 + lon_break2 = pgl_round(lon_break2 - 180); 5.728 + } else lon_dir2 = 0; 5.729 + /* iterate over all edges and vertices in cluster2 */ 5.730 + for (j2=0; j2<npoints2; j2++) { 5.731 + /* calculate index of next vertex */ 5.732 + k2 = (j2+1) % npoints2; 5.733 + /* skip last edge unless entry is (closed) outline or polygon */ 5.734 + if ( 5.735 + k2 == 0 && 5.736 + entrytype2 != PGL_ENTRY_OUTLINE && 5.737 + entrytype2 != PGL_ENTRY_POLYGON 5.738 + ) continue; 5.739 + /* use previously calculated values for lat1 and lon1 if possible */ 5.740 + if (j2) { 5.741 + lat21 = lat22; 5.742 + lon21 = lon22; 5.743 + } else { 5.744 + /* otherwise get latitude and longitude values of first vertex */ 5.745 + lat21 = points2[0].lat; 5.746 + lon21 = points2[0].lon; 5.747 + /* and consider longitude wrap-around for first vertex */ 5.748 + if (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360); 5.749 + else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360); 5.750 + } 5.751 + /* get latitude and longitude of next vertex */ 5.752 + lat22 = points2[k2].lat; 5.753 + lon22 = points2[k2].lon; 5.754 + /* consider longitude wrap-around for next vertex */ 5.755 + if (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360); 5.756 + else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360); 5.757 + /* skip degenerated edges */ 5.758 + if (lat21 == lat22 && lon21 == lon22) continue; 5.759 + /* perform another wrap-around where necessary */ 5.760 + /* TODO: improve performance of whole wrap-around mechanism */ 5.761 + wrapvalue = (lon21 + lon22) - (lon11 + lon12); 5.762 + if (wrapvalue > 360) { 5.763 + lon21 = pgl_round(lon21 - 360); 5.764 + lon22 = pgl_round(lon22 - 360); 5.765 + } else if (wrapvalue < -360) { 5.766 + lon21 = pgl_round(lon21 + 360); 5.767 + lon22 = pgl_round(lon22 + 360); 5.768 + } 5.769 + /* return true if segments overlap */ 5.770 + if ( 5.771 + pgl_lseg_crosses_line( 5.772 + lat11, lon11, lat12, lon12, 5.773 + lat21, lon21, lat22, lon22 5.774 + ) && pgl_lseg_crosses_line( 5.775 + lat21, lon21, lat22, lon22, 5.776 + lat11, lon11, lat12, lon12 5.777 + ) 5.778 + ) { 5.779 + return true; 5.780 + } 5.781 + } 5.782 + } 5.783 + } 5.784 + } 5.785 + /* otherwise return false */ 5.786 + return false; 5.787 +} 5.788 + 5.789 +/* check if second cluster is completely contained in first cluster */ 5.790 +static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) { 5.791 + if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false; 5.792 + if (pgl_any_cluster_points_in_cluster(inner, outer)) return false; 5.793 + if (pgl_outlines_overlap(outer, inner)) return false; 5.794 + return true; 5.795 +} 5.796 + 5.797 +/* check if two clusters overlap */ 5.798 +static bool pgl_clusters_overlap( 5.799 + pgl_cluster *cluster1, pgl_cluster *cluster2 5.800 +) { 5.801 + if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true; 5.802 + if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true; 5.803 + if (pgl_outlines_overlap(cluster1, cluster2)) return true; 5.804 + return false; 5.805 +} 5.806 + 5.807 + 5.808 +/* calculate (approximate) distance between point and cluster */ 5.809 +static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) { 5.810 + double comp; /* square of compression of meridians */ 5.811 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 5.812 + int entrytype; /* type of entry */ 5.813 + int npoints; /* number of points in entry */ 5.814 + pgl_point *points; /* array of points in entry */ 5.815 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 5.816 + double lon_break = 0; /* antipodal longitude of first vertex */ 5.817 + double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */ 5.818 + double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */ 5.819 + double lat0 = point->lat; /* latitude of point */ 5.820 + double lon0; /* (adjusted) longitude of point */ 5.821 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 5.822 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 5.823 + double s; /* scalar for vector calculations */ 5.824 + double dist; /* distance calculated in one step */ 5.825 + double min_dist = INFINITY; /* minimum distance */ 5.826 + /* distance is zero if point is contained in cluster */ 5.827 + if (pgl_point_in_cluster(point, cluster, false)) return 0; 5.828 + /* calculate approximate square compression of meridians */ 5.829 + comp = cos((lat0 / 180.0) * M_PI); 5.830 + comp *= comp; 5.831 + /* calculate exact square compression of meridians */ 5.832 + comp *= ( 5.833 + (1.0 - PGL_EPS2 * (1.0-comp)) * 5.834 + (1.0 - PGL_EPS2 * (1.0-comp)) / 5.835 + (PGL_SUBEPS2 * PGL_SUBEPS2) 5.836 + ); 5.837 + /* iterate over all entries */ 5.838 + for (i=0; i<cluster->nentries; i++) { 5.839 + /* get properties of entry */ 5.840 + entrytype = cluster->entries[i].entrytype; 5.841 + npoints = cluster->entries[i].npoints; 5.842 + points = PGL_ENTRY_POINTS(cluster, i); 5.843 + /* determine east/west orientation of first point of entry and calculate 5.844 + antipodal longitude */ 5.845 + lon_break = points[0].lon; 5.846 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 5.847 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 5.848 + else lon_dir = 0; 5.849 + /* determine covered longitude range */ 5.850 + for (j=0; j<npoints; j++) { 5.851 + /* get longitude of vertex */ 5.852 + lon1 = points[j].lon; 5.853 + /* adjust longitude to fix potential wrap-around */ 5.854 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 5.855 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 5.856 + /* update minimum and maximum longitude of polygon */ 5.857 + if (j == 0 || lon1 < lon_min) lon_min = lon1; 5.858 + if (j == 0 || lon1 > lon_max) lon_max = lon1; 5.859 + } 5.860 + /* adjust longitude wrap-around according to full longitude range */ 5.861 + lon_break = (lon_max + lon_min) / 2; 5.862 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 5.863 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 5.864 + /* get longitude of point */ 5.865 + lon0 = point->lon; 5.866 + /* consider longitude wrap-around for point */ 5.867 + if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360; 5.868 + else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360; 5.869 + /* iterate over all edges and vertices */ 5.870 + for (j=0; j<npoints; j++) { 5.871 + /* use previously calculated values for lat1 and lon1 if possible */ 5.872 + if (j) { 5.873 + lat1 = lat2; 5.874 + lon1 = lon2; 5.875 + } else { 5.876 + /* otherwise get latitude and longitude values of first vertex */ 5.877 + lat1 = points[0].lat; 5.878 + lon1 = points[0].lon; 5.879 + /* and consider longitude wrap-around for first vertex */ 5.880 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 5.881 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 5.882 + } 5.883 + /* calculate distance to vertex */ 5.884 + dist = pgl_distance(lat0, lon0, lat1, lon1); 5.885 + /* store calculated distance if smallest */ 5.886 + if (dist < min_dist) min_dist = dist; 5.887 + /* calculate index of next vertex */ 5.888 + k = (j+1) % npoints; 5.889 + /* skip last edge unless entry is (closed) outline or polygon */ 5.890 + if ( 5.891 + k == 0 && 5.892 + entrytype != PGL_ENTRY_OUTLINE && 5.893 + entrytype != PGL_ENTRY_POLYGON 5.894 + ) continue; 5.895 + /* get latitude and longitude of next vertex */ 5.896 + lat2 = points[k].lat; 5.897 + lon2 = points[k].lon; 5.898 + /* consider longitude wrap-around for next vertex */ 5.899 + if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360; 5.900 + else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360; 5.901 + /* go to next vertex and edge if edge is degenerated */ 5.902 + if (lat1 == lat2 && lon1 == lon2) continue; 5.903 + /* otherwise test if point can be projected onto edge of polygon */ 5.904 + s = ( 5.905 + ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) / 5.906 + ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1)) 5.907 + ); 5.908 + /* go to next vertex and edge if point cannot be projected */ 5.909 + if (!(s > 0 && s < 1)) continue; 5.910 + /* calculate distance from original point to projected point */ 5.911 + dist = pgl_distance( 5.912 + lat0, lon0, 5.913 + lat1 + s * (lat2-lat1), 5.914 + lon1 + s * (lon2-lon1) 5.915 + ); 5.916 + /* store calculated distance if smallest */ 5.917 + if (dist < min_dist) min_dist = dist; 5.918 + } 5.919 + } 5.920 + /* return minimum distance */ 5.921 + return min_dist; 5.922 +} 5.923 + 5.924 +/* calculate (approximate) distance between two clusters */ 5.925 +static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) { 5.926 + int i, j; /* i: entry, j: point in entry */ 5.927 + int npoints; /* number of points in entry */ 5.928 + pgl_point *points; /* array of points in entry */ 5.929 + double dist; /* distance calculated in one step */ 5.930 + double min_dist = INFINITY; /* minimum distance */ 5.931 + /* consider distance from each point in one cluster to the whole other */ 5.932 + for (i=0; i<cluster1->nentries; i++) { 5.933 + npoints = cluster1->entries[i].npoints; 5.934 + points = PGL_ENTRY_POINTS(cluster1, i); 5.935 + for (j=0; j<npoints; j++) { 5.936 + dist = pgl_point_cluster_distance(points+j, cluster2); 5.937 + if (dist == 0) return dist; 5.938 + if (dist < min_dist) min_dist = dist; 5.939 + } 5.940 + } 5.941 + /* consider distance from each point in other cluster to the first cluster */ 5.942 + for (i=0; i<cluster2->nentries; i++) { 5.943 + npoints = cluster2->entries[i].npoints; 5.944 + points = PGL_ENTRY_POINTS(cluster2, i); 5.945 + for (j=0; j<npoints; j++) { 5.946 + dist = pgl_point_cluster_distance(points+j, cluster1); 5.947 + if (dist == 0) return dist; 5.948 + if (dist < min_dist) min_dist = dist; 5.949 + } 5.950 + } 5.951 + return min_dist; 5.952 +} 5.953 + 5.954 +/* estimator function for distance between box and point */ 5.955 +/* always returns a smaller value than actually correct or zero */ 5.956 +static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) { 5.957 + double dlon; /* longitude range of box (delta longitude) */ 5.958 + double distance; /* return value */ 5.959 + /* return infinity if box is empty */ 5.960 + if (box->lat_min > box->lat_max) return INFINITY; 5.961 + /* return zero if point is inside box */ 5.962 + if (pgl_point_in_box(point, box)) return 0; 5.963 + /* calculate delta longitude */ 5.964 + dlon = box->lon_max - box->lon_min; 5.965 + if (dlon < 0) dlon += 360; /* 180th meridian crossed */ 5.966 + /* if delta longitude is greater than 150 degrees, perform safe fall-back */ 5.967 + if (dlon > 150) return 0; 5.968 + /* calculate lower limit for distance (formula below requires dlon <= 150) */ 5.969 + /* TODO: provide better estimation function to improve performance */ 5.970 + distance = ( 5.971 + (1.0-1e-14) * pgl_distance( 5.972 + point->lat, 5.973 + point->lon, 5.974 + (box->lat_min + box->lat_max) / 2, 5.975 + box->lon_min + dlon/2 5.976 + ) - pgl_distance( 5.977 + box->lat_min, box->lon_min, 5.978 + box->lat_max, box->lon_max 5.979 + ) 5.980 + ); 5.981 + /* truncate negative results to zero */ 5.982 + if (distance <= 0) distance = 0; 5.983 + /* return result */ 5.984 + return distance; 5.985 +} 5.986 + 5.987 + 5.988 +/*-------------------------------* 5.989 + * Monte Carlo based functions * 5.990 + *-------------------------------*/ 5.991 + 5.992 +/* half of (spherical) earth's surface area */ 5.993 +#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI) 5.994 + 5.995 +/* golden angle in radians */ 5.996 +#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0)) 5.997 + 5.998 +/* create a list of sample points covering a bounding circle 5.999 + and return covered area */ 5.1000 +static double pgl_sample_points( 5.1001 + pgl_point *center, /* center of bounding circle */ 5.1002 + double radius, /* radius of bounding circle */ 5.1003 + int samples, /* number of sample points */ 5.1004 + pgl_point *result /* pointer to result array */ 5.1005 +) { 5.1006 + double double_share = 2.0; /* double of covered share of earth's surface */ 5.1007 + double double_share_div_samples; /* double_share divided by sample count */ 5.1008 + int i; 5.1009 + double t; /* parameter of spiral laid on (spherical) earth's surface */ 5.1010 + double x, y, z; /* normalized coordinates of point on non-rotated spiral */ 5.1011 + double sin_phi; /* sine of sph. coordinate of point of non-rotated spiral */ 5.1012 + double lambda; /* other sph. coordinate of point of non-rotated spiral */ 5.1013 + double rot = (0.5 - center->lat / 180.0) * M_PI; /* needed rot. (in rad) */ 5.1014 + double cos_rot = cos(rot); /* cosine of rotation by latitude */ 5.1015 + double sin_rot = sin(rot); /* sine of rotation by latitude */ 5.1016 + double x_rot, z_rot; /* normalized coordinates of point on rotated spiral */ 5.1017 + double center_lon = center->lon; /* second rotation in degree */ 5.1018 + /* add safety margin to bounding circle because of spherical approximation */ 5.1019 + radius *= PGL_SPHEROID_A / PGL_RADIUS; 5.1020 + /* if whole earth is covered, use initialized value, otherwise calculate 5.1021 + share of covered area (multiplied by 2) */ 5.1022 + if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS); 5.1023 + /* divide double_share by sample count for later calculations */ 5.1024 + double_share_div_samples = double_share / samples; 5.1025 + /* generate sample points */ 5.1026 + for (i=0; i<samples; i++) { 5.1027 + /* use an offset of 1/2 to avoid too dense clustering at spiral center */ 5.1028 + t = 0.5 + i; 5.1029 + /* calculate normalized coordinates of point on non-rotated spiral */ 5.1030 + z = 1.0 - double_share_div_samples * t; 5.1031 + sin_phi = sqrt(1.0 - z*z); 5.1032 + lambda = t * PGL_GOLDEN_ANGLE; 5.1033 + x = sin_phi * cos(lambda); 5.1034 + y = sin_phi * sin(lambda); 5.1035 + /* rotate spiral by latitude value of bounding circle */ 5.1036 + x_rot = cos_rot * x + sin_rot * z; 5.1037 + z_rot = cos_rot * z - sin_rot * x; 5.1038 + /* set resulting sample point in result array */ 5.1039 + /* (while performing second rotation by bounding circle longitude) */ 5.1040 + result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI); 5.1041 + result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI); 5.1042 + } 5.1043 + /* return covered area */ 5.1044 + return PGL_HALF_SURFACE * double_share; 5.1045 +} 5.1046 + 5.1047 +/* calculate covered area using Monte Carlo method */ 5.1048 +/* TODO: inefficient, should be replaced by different method */ 5.1049 +static double pgl_monte_carlo_area(pgl_cluster *cluster, int samples) { 5.1050 + pgl_point *points; 5.1051 + double area; 5.1052 + int i; 5.1053 + int matches = 0; 5.1054 + if (samples > PGL_CLUSTER_MAXPOINTS) { 5.1055 + ereport(ERROR, ( 5.1056 + errcode(ERRCODE_DATA_EXCEPTION), 5.1057 + errmsg( 5.1058 + "too many sample points for Monte Carlo method (maximum %i)", 5.1059 + PGL_CLUSTER_MAXPOINTS 5.1060 + ) 5.1061 + )); 5.1062 + } 5.1063 + if (!(cluster->bounding.radius > 0)) return 0; 5.1064 + for (i=0; ; i++) { 5.1065 + if (i == cluster->nentries) return 0; 5.1066 + if (cluster->entries[i].entrytype == PGL_ENTRY_POLYGON) break; 5.1067 + } 5.1068 + points = palloc(samples * sizeof(pgl_point)); 5.1069 + area = pgl_sample_points( 5.1070 + &cluster->bounding.center, 5.1071 + cluster->bounding.radius, 5.1072 + samples, 5.1073 + points 5.1074 + ); 5.1075 + for (i=0; i<samples; i++) { 5.1076 + if (pgl_point_in_cluster(points+i, cluster, true)) matches++; 5.1077 + } 5.1078 + pfree(points); 5.1079 + return area * ((double)matches / samples); 5.1080 +} 5.1081 + 5.1082 +/* fair distance between point and cluster (see README file for explanation) */ 5.1083 +static double pgl_fair_distance( 5.1084 + pgl_point *point, pgl_cluster *cluster, int samples 5.1085 +) { 5.1086 + double distance; /* shortest distance from point to cluster */ 5.1087 + pgl_point *points; /* sample points for Monte Carlo method */ 5.1088 + double area; /* area covered by sample points */ 5.1089 + int i; 5.1090 + int matches = 0; /* number of points within enlarged area (multiplied 5.1091 + by two to be able to store half matches) */ 5.1092 + double result; /* result determined by Monte Carlo method */ 5.1093 + /* limit sample count to avoid integer overflows on memory allocation */ 5.1094 + if (samples > PGL_CLUSTER_MAXPOINTS) { 5.1095 + ereport(ERROR, ( 5.1096 + errcode(ERRCODE_DATA_EXCEPTION), 5.1097 + errmsg( 5.1098 + "too many sample points for Monte Carlo method (maximum %i)", 5.1099 + PGL_CLUSTER_MAXPOINTS 5.1100 + ) 5.1101 + )); 5.1102 + } 5.1103 + /* calculate shortest distance from point to cluster */ 5.1104 + distance = pgl_point_cluster_distance(point, cluster); 5.1105 + /* if cluster consists of a single point or has no bounding circle with 5.1106 + positive radius, simply return distance */ 5.1107 + if ( 5.1108 + (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) || 5.1109 + !(cluster->bounding.radius > 0) 5.1110 + ) return distance; 5.1111 + /* if cluster consists of two points which are twice as far apart, return 5.1112 + distance between point and cluster multiplied by square root of two */ 5.1113 + if ( 5.1114 + cluster->nentries == 2 && 5.1115 + cluster->entries[0].entrytype == PGL_ENTRY_POINT && 5.1116 + cluster->entries[1].entrytype == PGL_ENTRY_POINT && 5.1117 + pgl_distance( 5.1118 + PGL_ENTRY_POINTS(cluster, 0)[0].lat, 5.1119 + PGL_ENTRY_POINTS(cluster, 0)[0].lon, 5.1120 + PGL_ENTRY_POINTS(cluster, 1)[0].lat, 5.1121 + PGL_ENTRY_POINTS(cluster, 1)[0].lon 5.1122 + ) >= 2.0 * distance 5.1123 + ) { 5.1124 + return distance * M_SQRT2; 5.1125 + } 5.1126 + /* otherwise create sample points for Monte Carlo method and determine area 5.1127 + covered by sample points */ 5.1128 + points = palloc(samples * sizeof(pgl_point)); 5.1129 + area = pgl_sample_points( 5.1130 + &cluster->bounding.center, 5.1131 + cluster->bounding.radius + distance, /* pad bounding circle by distance */ 5.1132 + samples, 5.1133 + points 5.1134 + ); 5.1135 + /* perform Monte Carlo method */ 5.1136 + if (distance > 0) { 5.1137 + /* point is outside cluster */ 5.1138 + for (i=0; i<samples; i++) { 5.1139 + /* count sample poitns within cluster as half match */ 5.1140 + if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1; 5.1141 + /* count sample points outside of cluster but within cluster grown by 5.1142 + distance as full match */ 5.1143 + else if ( 5.1144 + pgl_point_cluster_distance(points+i, cluster) < distance 5.1145 + ) matches += 2; 5.1146 + } 5.1147 + } else { 5.1148 + /* if point is within cluster, 5.1149 + just count sample points within cluster as half match */ 5.1150 + for (i=0; i<samples; i++) { 5.1151 + if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1; 5.1152 + } 5.1153 + } 5.1154 + /* release memory for sample points needed by Monte Carlo method */ 5.1155 + pfree(points); 5.1156 + /* convert area determined by Monte Carlo method into distance 5.1157 + (i.e. radius of circle with the same area) */ 5.1158 + result = sqrt(area * ((double)matches / 2.0 / samples) / M_PI); 5.1159 + /* return result only if it is greater than the distance between point and 5.1160 + cluster to avoid unexpected results because of errors due to limited 5.1161 + precision */ 5.1162 + if (result > distance) return result; 5.1163 + /* otherwise return distance between point and cluster */ 5.1164 + else return distance; 5.1165 +} 5.1166 + 5.1167 + 5.1168 +/*-------------------------------------------------* 5.1169 + * geographic index based on space-filling curve * 5.1170 + *-------------------------------------------------*/ 5.1171 + 5.1172 +/* number of bytes used for geographic (center) position in keys */ 5.1173 +#define PGL_KEY_LATLON_BYTELEN 7 5.1174 + 5.1175 +/* maximum reference value for logarithmic size of geographic objects */ 5.1176 +#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */ 5.1177 + 5.1178 +/* pointer to index key (either pgl_pointkey or pgl_areakey) */ 5.1179 +typedef unsigned char *pgl_keyptr; 5.1180 + 5.1181 +/* index key for points (objects with zero area) on the spheroid */ 5.1182 +/* bit 0..55: interspersed bits of latitude and longitude, 5.1183 + bit 56..57: always zero, 5.1184 + bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */ 5.1185 +typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1]; 5.1186 + 5.1187 +/* index key for geographic objects on spheroid with area greater than zero */ 5.1188 +/* bit 0..55: interspersed bits of latitude and longitude of center point, 5.1189 + bit 56: always set to 1, 5.1190 + bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.), 5.1191 + bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to 5.1192 + PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth 5.1193 + = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL 5.1194 + (with interspersed bits = 0 and node depth = 0) for keys which 5.1195 + cover both empty and non-empty objects */ 5.1196 + 5.1197 +typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2]; 5.1198 + 5.1199 +/* helper macros for reading/writing index keys */ 5.1200 +#define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN 5.1201 +#define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1) 5.1202 +#define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8) 5.1203 +#define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1) 5.1204 +#define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1) 5.1205 +#define PGL_AREAKEY_TYPEMASK 0x80 5.1206 +#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8))) 5.1207 +#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \ 5.1208 + ( PGL_KEY_LATLONBIT(key1, n) ^ \ 5.1209 + PGL_KEY_LATLONBIT(key2, n) ) 5.1210 +#define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 5.1211 + PGL_AREAKEY_TYPEMASK) 5.1212 +#define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 5.1213 + (PGL_AREAKEY_TYPEMASK-1)) 5.1214 +#define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET]) 5.1215 +#define PGL_KEY_OBJSIZE_EMPTY 126 5.1216 +#define PGL_KEY_OBJSIZE_UNIVERSAL 127 5.1217 +#define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \ 5.1218 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 5.1219 + PGL_KEY_OBJSIZE_EMPTY ) 5.1220 +#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \ 5.1221 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 5.1222 + PGL_KEY_OBJSIZE_UNIVERSAL ) 5.1223 + 5.1224 +/* set area key to match empty objects only */ 5.1225 +static void pgl_key_set_empty(pgl_keyptr key) { 5.1226 + memset(key, 0, sizeof(pgl_areakey)); 5.1227 + /* Note: setting node depth to maximum is required for picksplit function */ 5.1228 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 5.1229 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY; 5.1230 +} 5.1231 + 5.1232 +/* set area key to match any object (including empty objects) */ 5.1233 +static void pgl_key_set_universal(pgl_keyptr key) { 5.1234 + memset(key, 0, sizeof(pgl_areakey)); 5.1235 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK; 5.1236 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL; 5.1237 +} 5.1238 + 5.1239 +/* convert a point on earth into a max-depth key to be used in index */ 5.1240 +static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) { 5.1241 + double lat = point->lat; 5.1242 + double lon = point->lon; 5.1243 + int i; 5.1244 + /* clear latitude and longitude bits */ 5.1245 + memset(key, 0, PGL_KEY_LATLON_BYTELEN); 5.1246 + /* set node depth to maximum and type bit to zero */ 5.1247 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH; 5.1248 + /* iterate over all latitude/longitude bit pairs */ 5.1249 + for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) { 5.1250 + /* determine latitude bit */ 5.1251 + if (lat >= 0) { 5.1252 + key[i/4] |= 0x80 >> (2*(i%4)); 5.1253 + lat *= 2; lat -= 90; 5.1254 + } else { 5.1255 + lat *= 2; lat += 90; 5.1256 + } 5.1257 + /* determine longitude bit */ 5.1258 + if (lon >= 0) { 5.1259 + key[i/4] |= 0x80 >> (2*(i%4)+1); 5.1260 + lon *= 2; lon -= 180; 5.1261 + } else { 5.1262 + lon *= 2; lon += 180; 5.1263 + } 5.1264 + } 5.1265 +} 5.1266 + 5.1267 +/* convert a circle on earth into a max-depth key to be used in an index */ 5.1268 +static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) { 5.1269 + /* handle special case of empty circle */ 5.1270 + if (circle->radius < 0) { 5.1271 + pgl_key_set_empty(key); 5.1272 + return; 5.1273 + } 5.1274 + /* perform same action as for point keys */ 5.1275 + pgl_point_to_key(&(circle->center), key); 5.1276 + /* but overwrite type and node depth to fit area index key */ 5.1277 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 5.1278 + /* check if radius is greater than (or equal to) reference size */ 5.1279 + /* (treat equal values as greater values for numerical safety) */ 5.1280 + if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) { 5.1281 + /* if yes, set logarithmic size to zero */ 5.1282 + key[PGL_KEY_OBJSIZE_OFFSET] = 0; 5.1283 + } else { 5.1284 + /* otherwise, determine logarithmic size iteratively */ 5.1285 + /* (one step is equivalent to a factor of sqrt(2)) */ 5.1286 + double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2; 5.1287 + int objsize = 1; 5.1288 + while (objsize < PGL_AREAKEY_MAXOBJSIZE) { 5.1289 + /* stop when radius is greater than (or equal to) adjusted reference */ 5.1290 + /* (treat equal values as greater values for numerical safety) */ 5.1291 + if (circle->radius >= reference) break; 5.1292 + reference /= M_SQRT2; 5.1293 + objsize++; 5.1294 + } 5.1295 + /* set logarithmic size to determined value */ 5.1296 + key[PGL_KEY_OBJSIZE_OFFSET] = objsize; 5.1297 + } 5.1298 +} 5.1299 + 5.1300 +/* check if one key is subkey of another key or vice versa */ 5.1301 +static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) { 5.1302 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 5.1303 + /* determine smallest depth */ 5.1304 + int depth1 = PGL_KEY_NODEDEPTH(key1); 5.1305 + int depth2 = PGL_KEY_NODEDEPTH(key2); 5.1306 + int depth = (depth1 < depth2) ? depth1 : depth2; 5.1307 + /* check if keys are area keys (assuming that both keys have same type) */ 5.1308 + if (PGL_KEY_IS_AREAKEY(key1)) { 5.1309 + int j = 0; /* bit offset for logarithmic object size bits */ 5.1310 + int k = 0; /* bit offset for latitude and longitude */ 5.1311 + /* fetch logarithmic object size information */ 5.1312 + int objsize1 = PGL_KEY_OBJSIZE(key1); 5.1313 + int objsize2 = PGL_KEY_OBJSIZE(key2); 5.1314 + /* handle special cases for empty objects (universal and empty keys) */ 5.1315 + if ( 5.1316 + objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL || 5.1317 + objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL 5.1318 + ) return true; 5.1319 + if ( 5.1320 + objsize1 == PGL_KEY_OBJSIZE_EMPTY || 5.1321 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 5.1322 + ) return objsize1 == objsize2; 5.1323 + /* iterate through key bits */ 5.1324 + for (i=0; i<depth; i++) { 5.1325 + /* every second bit is a bit describing the object size */ 5.1326 + if (i%2 == 0) { 5.1327 + /* check if object size bit is different in both keys (objsize1 and 5.1328 + objsize2 describe the minimum index when object size bit is set) */ 5.1329 + if ( 5.1330 + (objsize1 <= j && objsize2 > j) || 5.1331 + (objsize2 <= j && objsize1 > j) 5.1332 + ) { 5.1333 + /* bit differs, therefore keys are in separate branches */ 5.1334 + return false; 5.1335 + } 5.1336 + /* increase bit counter for object size bits */ 5.1337 + j++; 5.1338 + } 5.1339 + /* all other bits describe latitude and longitude */ 5.1340 + else { 5.1341 + /* check if bit differs in both keys */ 5.1342 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) { 5.1343 + /* bit differs, therefore keys are in separate branches */ 5.1344 + return false; 5.1345 + } 5.1346 + /* increase bit counter for latitude/longitude bits */ 5.1347 + k++; 5.1348 + } 5.1349 + } 5.1350 + } 5.1351 + /* if not, keys are point keys */ 5.1352 + else { 5.1353 + /* iterate through key bits */ 5.1354 + for (i=0; i<depth; i++) { 5.1355 + /* check if bit differs in both keys */ 5.1356 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) { 5.1357 + /* bit differs, therefore keys are in separate branches */ 5.1358 + return false; 5.1359 + } 5.1360 + } 5.1361 + } 5.1362 + /* return true because keys are in the same branch */ 5.1363 + return true; 5.1364 +} 5.1365 + 5.1366 +/* combine two keys into new key which covers both original keys */ 5.1367 +/* (result stored in first argument) */ 5.1368 +static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) { 5.1369 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 5.1370 + /* determine smallest depth */ 5.1371 + int depth1 = PGL_KEY_NODEDEPTH(dst); 5.1372 + int depth2 = PGL_KEY_NODEDEPTH(src); 5.1373 + int depth = (depth1 < depth2) ? depth1 : depth2; 5.1374 + /* check if keys are area keys (assuming that both keys have same type) */ 5.1375 + if (PGL_KEY_IS_AREAKEY(dst)) { 5.1376 + pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */ 5.1377 + int j = 0; /* bit offset for logarithmic object size bits */ 5.1378 + int k = 0; /* bit offset for latitude and longitude */ 5.1379 + /* fetch logarithmic object size information */ 5.1380 + int objsize1 = PGL_KEY_OBJSIZE(dst); 5.1381 + int objsize2 = PGL_KEY_OBJSIZE(src); 5.1382 + /* handle special cases for empty objects (universal and empty keys) */ 5.1383 + if ( 5.1384 + objsize1 > PGL_AREAKEY_MAXOBJSIZE || 5.1385 + objsize2 > PGL_AREAKEY_MAXOBJSIZE 5.1386 + ) { 5.1387 + if ( 5.1388 + objsize1 == PGL_KEY_OBJSIZE_EMPTY && 5.1389 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 5.1390 + ) pgl_key_set_empty(dst); 5.1391 + else pgl_key_set_universal(dst); 5.1392 + return; 5.1393 + } 5.1394 + /* iterate through key bits */ 5.1395 + for (i=0; i<depth; i++) { 5.1396 + /* every second bit is a bit describing the object size */ 5.1397 + if (i%2 == 0) { 5.1398 + /* increase bit counter for object size bits first */ 5.1399 + /* (handy when setting objsize variable) */ 5.1400 + j++; 5.1401 + /* check if object size bit is set in neither key */ 5.1402 + if (objsize1 >= j && objsize2 >= j) { 5.1403 + /* set objsize in destination buffer to indicate that size bit is 5.1404 + unset in destination buffer at the current bit position */ 5.1405 + dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j; 5.1406 + } 5.1407 + /* break if object size bit is set in one key only */ 5.1408 + else if (objsize1 >= j || objsize2 >= j) break; 5.1409 + } 5.1410 + /* all other bits describe latitude and longitude */ 5.1411 + else { 5.1412 + /* break if bit differs in both keys */ 5.1413 + if (PGL_KEY_LATLONBIT(dst, k)) { 5.1414 + if (!PGL_KEY_LATLONBIT(src, k)) break; 5.1415 + /* but set bit in destination buffer if bit is set in both keys */ 5.1416 + dstbuf[k/8] |= 0x80 >> (k%8); 5.1417 + } else if (PGL_KEY_LATLONBIT(src, k)) break; 5.1418 + /* increase bit counter for latitude/longitude bits */ 5.1419 + k++; 5.1420 + } 5.1421 + } 5.1422 + /* set common node depth and type bit (type bit = 1) */ 5.1423 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i; 5.1424 + /* copy contents of destination buffer to first key */ 5.1425 + memcpy(dst, dstbuf, sizeof(pgl_areakey)); 5.1426 + } 5.1427 + /* if not, keys are point keys */ 5.1428 + else { 5.1429 + pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */ 5.1430 + /* iterate through key bits */ 5.1431 + for (i=0; i<depth; i++) { 5.1432 + /* break if bit differs in both keys */ 5.1433 + if (PGL_KEY_LATLONBIT(dst, i)) { 5.1434 + if (!PGL_KEY_LATLONBIT(src, i)) break; 5.1435 + /* but set bit in destination buffer if bit is set in both keys */ 5.1436 + dstbuf[i/8] |= 0x80 >> (i%8); 5.1437 + } else if (PGL_KEY_LATLONBIT(src, i)) break; 5.1438 + } 5.1439 + /* set common node depth (type bit = 0) */ 5.1440 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i; 5.1441 + /* copy contents of destination buffer to first key */ 5.1442 + memcpy(dst, dstbuf, sizeof(pgl_pointkey)); 5.1443 + } 5.1444 +} 5.1445 + 5.1446 +/* determine center(!) boundaries and radius estimation of index key */ 5.1447 +static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) { 5.1448 + int i; 5.1449 + /* determine node depth */ 5.1450 + int depth = PGL_KEY_NODEDEPTH(key); 5.1451 + /* center point of possible result */ 5.1452 + double lat = 0; 5.1453 + double lon = 0; 5.1454 + /* maximum distance of real center point from key center */ 5.1455 + double dlat = 90; 5.1456 + double dlon = 180; 5.1457 + /* maximum radius of contained objects */ 5.1458 + double radius = 0; /* always return zero for point index keys */ 5.1459 + /* check if key is area key */ 5.1460 + if (PGL_KEY_IS_AREAKEY(key)) { 5.1461 + /* get logarithmic object size */ 5.1462 + int objsize = PGL_KEY_OBJSIZE(key); 5.1463 + /* handle special cases for empty objects (universal and empty keys) */ 5.1464 + if (objsize == PGL_KEY_OBJSIZE_EMPTY) { 5.1465 + pgl_box_set_empty(box); 5.1466 + return 0; 5.1467 + } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) { 5.1468 + box->lat_min = -90; 5.1469 + box->lat_max = 90; 5.1470 + box->lon_min = -180; 5.1471 + box->lon_max = 180; 5.1472 + return 0; /* any value >= 0 would do */ 5.1473 + } 5.1474 + /* calculate maximum possible radius of objects covered by the given key */ 5.1475 + if (objsize == 0) radius = INFINITY; 5.1476 + else { 5.1477 + radius = PGL_AREAKEY_REFOBJSIZE; 5.1478 + while (--objsize) radius /= M_SQRT2; 5.1479 + } 5.1480 + /* iterate over latitude and longitude bits in key */ 5.1481 + /* (every second bit is a latitude or longitude bit) */ 5.1482 + for (i=0; i<depth/2; i++) { 5.1483 + /* check if latitude bit */ 5.1484 + if (i%2 == 0) { 5.1485 + /* cut latitude dimension in half */ 5.1486 + dlat /= 2; 5.1487 + /* increase center latitude if bit is 1, otherwise decrease */ 5.1488 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 5.1489 + else lat -= dlat; 5.1490 + } 5.1491 + /* otherwise longitude bit */ 5.1492 + else { 5.1493 + /* cut longitude dimension in half */ 5.1494 + dlon /= 2; 5.1495 + /* increase center longitude if bit is 1, otherwise decrease */ 5.1496 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 5.1497 + else lon -= dlon; 5.1498 + } 5.1499 + } 5.1500 + } 5.1501 + /* if not, keys are point keys */ 5.1502 + else { 5.1503 + /* iterate over all bits in key */ 5.1504 + for (i=0; i<depth; i++) { 5.1505 + /* check if latitude bit */ 5.1506 + if (i%2 == 0) { 5.1507 + /* cut latitude dimension in half */ 5.1508 + dlat /= 2; 5.1509 + /* increase center latitude if bit is 1, otherwise decrease */ 5.1510 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 5.1511 + else lat -= dlat; 5.1512 + } 5.1513 + /* otherwise longitude bit */ 5.1514 + else { 5.1515 + /* cut longitude dimension in half */ 5.1516 + dlon /= 2; 5.1517 + /* increase center longitude if bit is 1, otherwise decrease */ 5.1518 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 5.1519 + else lon -= dlon; 5.1520 + } 5.1521 + } 5.1522 + } 5.1523 + /* calculate boundaries from center point and remaining dlat and dlon */ 5.1524 + /* (return values through pointer to box) */ 5.1525 + box->lat_min = lat - dlat; 5.1526 + box->lat_max = lat + dlat; 5.1527 + box->lon_min = lon - dlon; 5.1528 + box->lon_max = lon + dlon; 5.1529 + /* return radius (as a function return value) */ 5.1530 + return radius; 5.1531 +} 5.1532 + 5.1533 +/* estimator function for distance between point and index key */ 5.1534 +/* always returns a smaller value than actually correct or zero */ 5.1535 +static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) { 5.1536 + pgl_box box; /* center(!) bounding box of area index key */ 5.1537 + /* calculate center(!) bounding box and maximum radius of objects covered 5.1538 + by area index key (radius is zero for point index keys) */ 5.1539 + double distance = pgl_key_to_box(key, &box); 5.1540 + /* calculate estimated distance between bounding box of center point of 5.1541 + indexed object and point passed as second argument, then substract maximum 5.1542 + radius of objects covered by index key */ 5.1543 + distance = pgl_estimate_point_box_distance(point, &box) - distance; 5.1544 + /* truncate negative results to zero */ 5.1545 + if (distance <= 0) distance = 0; 5.1546 + /* return result */ 5.1547 + return distance; 5.1548 +} 5.1549 + 5.1550 + 5.1551 +/*---------------------------------* 5.1552 + * helper functions for text I/O * 5.1553 + *---------------------------------*/ 5.1554 + 5.1555 +#define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */ 5.1556 + 5.1557 +/* convert floating point number to string (round-trip safe) */ 5.1558 +static void pgl_print_float(char *buf, double flt) { 5.1559 + /* check if number is integral */ 5.1560 + if (trunc(flt) == flt) { 5.1561 + /* for integral floats use maximum precision */ 5.1562 + snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 5.1563 + } else { 5.1564 + /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */ 5.1565 + snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt); 5.1566 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt); 5.1567 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 5.1568 + } 5.1569 +} 5.1570 + 5.1571 +/* convert latitude floating point number (in degrees) to string */ 5.1572 +static void pgl_print_lat(char *buf, double lat) { 5.1573 + if (signbit(lat)) { 5.1574 + /* treat negative latitudes (including -0) as south */ 5.1575 + snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat); 5.1576 + } else { 5.1577 + /* treat positive latitudes (including +0) as north */ 5.1578 + snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat); 5.1579 + } 5.1580 +} 5.1581 + 5.1582 +/* convert longitude floating point number (in degrees) to string */ 5.1583 +static void pgl_print_lon(char *buf, double lon) { 5.1584 + if (signbit(lon)) { 5.1585 + /* treat negative longitudes (including -0) as west */ 5.1586 + snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon); 5.1587 + } else { 5.1588 + /* treat positive longitudes (including +0) as east */ 5.1589 + snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon); 5.1590 + } 5.1591 +} 5.1592 + 5.1593 +/* bit masks used as return value of pgl_scan() function */ 5.1594 +#define PGL_SCAN_NONE 0 /* no value has been parsed */ 5.1595 +#define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */ 5.1596 +#define PGL_SCAN_LON (1<<1) /* longitude has been parsed */ 5.1597 +#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */ 5.1598 + 5.1599 +/* parse a coordinate (can be latitude or longitude) */ 5.1600 +static int pgl_scan(char **str, double *lat, double *lon) { 5.1601 + double val; 5.1602 + int len; 5.1603 + if ( 5.1604 + sscanf(*str, " N %lf %n", &val, &len) || 5.1605 + sscanf(*str, " n %lf %n", &val, &len) 5.1606 + ) { 5.1607 + *str += len; *lat = val; return PGL_SCAN_LAT; 5.1608 + } 5.1609 + if ( 5.1610 + sscanf(*str, " S %lf %n", &val, &len) || 5.1611 + sscanf(*str, " s %lf %n", &val, &len) 5.1612 + ) { 5.1613 + *str += len; *lat = -val; return PGL_SCAN_LAT; 5.1614 + } 5.1615 + if ( 5.1616 + sscanf(*str, " E %lf %n", &val, &len) || 5.1617 + sscanf(*str, " e %lf %n", &val, &len) 5.1618 + ) { 5.1619 + *str += len; *lon = val; return PGL_SCAN_LON; 5.1620 + } 5.1621 + if ( 5.1622 + sscanf(*str, " W %lf %n", &val, &len) || 5.1623 + sscanf(*str, " w %lf %n", &val, &len) 5.1624 + ) { 5.1625 + *str += len; *lon = -val; return PGL_SCAN_LON; 5.1626 + } 5.1627 + return PGL_SCAN_NONE; 5.1628 +} 5.1629 + 5.1630 + 5.1631 +/*-----------------* 5.1632 + * SQL functions * 5.1633 + *-----------------*/ 5.1634 + 5.1635 +/* Note: These function names use "epoint", "ebox", etc. notation here instead 5.1636 + of "point", "box", etc. in order to distinguish them from any previously 5.1637 + defined functions. */ 5.1638 + 5.1639 +/* function needed for dummy types and/or not implemented features */ 5.1640 +PG_FUNCTION_INFO_V1(pgl_notimpl); 5.1641 +Datum pgl_notimpl(PG_FUNCTION_ARGS) { 5.1642 + ereport(ERROR, (errmsg("not implemented by pgLatLon"))); 5.1643 +} 5.1644 + 5.1645 +/* set point to latitude and longitude (including checks) */ 5.1646 +static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) { 5.1647 + /* reject infinite or NaN values */ 5.1648 + if (!isfinite(lat) || !isfinite(lon)) { 5.1649 + ereport(ERROR, ( 5.1650 + errcode(ERRCODE_DATA_EXCEPTION), 5.1651 + errmsg("epoint requires finite coordinates") 5.1652 + )); 5.1653 + } 5.1654 + /* check latitude bounds */ 5.1655 + if (lat < -90) { 5.1656 + ereport(WARNING, (errmsg("latitude exceeds south pole"))); 5.1657 + lat = -90; 5.1658 + } else if (lat > 90) { 5.1659 + ereport(WARNING, (errmsg("latitude exceeds north pole"))); 5.1660 + lat = 90; 5.1661 + } 5.1662 + /* check longitude bounds */ 5.1663 + if (lon < -180) { 5.1664 + ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized"))); 5.1665 + lon += 360 - trunc(lon / 360) * 360; 5.1666 + } else if (lon > 180) { 5.1667 + ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized"))); 5.1668 + lon -= 360 + trunc(lon / 360) * 360; 5.1669 + } 5.1670 + /* store rounded latitude/longitude values for round-trip safety */ 5.1671 + point->lat = pgl_round(lat); 5.1672 + point->lon = pgl_round(lon); 5.1673 +} 5.1674 + 5.1675 +/* create point ("epoint" in SQL) from latitude and longitude */ 5.1676 +PG_FUNCTION_INFO_V1(pgl_create_epoint); 5.1677 +Datum pgl_create_epoint(PG_FUNCTION_ARGS) { 5.1678 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 5.1679 + pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1)); 5.1680 + PG_RETURN_POINTER(point); 5.1681 +} 5.1682 + 5.1683 +/* parse point ("epoint" in SQL) */ 5.1684 +/* format: '[NS]<float> [EW]<float>' */ 5.1685 +PG_FUNCTION_INFO_V1(pgl_epoint_in); 5.1686 +Datum pgl_epoint_in(PG_FUNCTION_ARGS) { 5.1687 + char *str = PG_GETARG_CSTRING(0); /* input string */ 5.1688 + char *strptr = str; /* current position within string */ 5.1689 + int done = 0; /* bit mask storing if latitude or longitude was read */ 5.1690 + double lat, lon; /* parsed values as double precision floats */ 5.1691 + pgl_point *point; /* return value (to be palloc'ed) */ 5.1692 + /* parse two floats (each latitude or longitude) separated by white-space */ 5.1693 + done |= pgl_scan(&strptr, &lat, &lon); 5.1694 + if (strptr != str && isspace(strptr[-1])) { 5.1695 + done |= pgl_scan(&strptr, &lat, &lon); 5.1696 + } 5.1697 + /* require end of string, and latitude and longitude parsed successfully */ 5.1698 + if (strptr[0] || done != PGL_SCAN_LATLON) { 5.1699 + ereport(ERROR, ( 5.1700 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.1701 + errmsg("invalid input syntax for type epoint: \"%s\"", str) 5.1702 + )); 5.1703 + } 5.1704 + /* allocate memory for result */ 5.1705 + point = (pgl_point *)palloc(sizeof(pgl_point)); 5.1706 + /* set latitude and longitude (and perform checks) */ 5.1707 + pgl_epoint_set_latlon(point, lat, lon); 5.1708 + /* return result */ 5.1709 + PG_RETURN_POINTER(point); 5.1710 +} 5.1711 + 5.1712 +/* create box ("ebox" in SQL) that is empty */ 5.1713 +PG_FUNCTION_INFO_V1(pgl_create_empty_ebox); 5.1714 +Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) { 5.1715 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 5.1716 + pgl_box_set_empty(box); 5.1717 + PG_RETURN_POINTER(box); 5.1718 +} 5.1719 + 5.1720 +/* set box to given boundaries (including checks) */ 5.1721 +static void pgl_ebox_set_boundaries( 5.1722 + pgl_box *box, 5.1723 + double lat_min, double lat_max, double lon_min, double lon_max 5.1724 +) { 5.1725 + /* if minimum latitude is greater than maximum latitude, return empty box */ 5.1726 + if (lat_min > lat_max) { 5.1727 + pgl_box_set_empty(box); 5.1728 + return; 5.1729 + } 5.1730 + /* otherwise reject infinite or NaN values */ 5.1731 + if ( 5.1732 + !isfinite(lat_min) || !isfinite(lat_max) || 5.1733 + !isfinite(lon_min) || !isfinite(lon_max) 5.1734 + ) { 5.1735 + ereport(ERROR, ( 5.1736 + errcode(ERRCODE_DATA_EXCEPTION), 5.1737 + errmsg("ebox requires finite coordinates") 5.1738 + )); 5.1739 + } 5.1740 + /* check latitude bounds */ 5.1741 + if (lat_max < -90) { 5.1742 + ereport(WARNING, (errmsg("northern latitude exceeds south pole"))); 5.1743 + lat_max = -90; 5.1744 + } else if (lat_max > 90) { 5.1745 + ereport(WARNING, (errmsg("northern latitude exceeds north pole"))); 5.1746 + lat_max = 90; 5.1747 + } 5.1748 + if (lat_min < -90) { 5.1749 + ereport(WARNING, (errmsg("southern latitude exceeds south pole"))); 5.1750 + lat_min = -90; 5.1751 + } else if (lat_min > 90) { 5.1752 + ereport(WARNING, (errmsg("southern latitude exceeds north pole"))); 5.1753 + lat_min = 90; 5.1754 + } 5.1755 + /* check if all longitudes are included */ 5.1756 + if (lon_max - lon_min >= 360) { 5.1757 + if (lon_max - lon_min > 360) ereport(WARNING, ( 5.1758 + errmsg("longitude coverage greater than 360 degrees") 5.1759 + )); 5.1760 + lon_min = -180; 5.1761 + lon_max = 180; 5.1762 + } else { 5.1763 + /* normalize longitude bounds */ 5.1764 + if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360; 5.1765 + else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360; 5.1766 + if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360; 5.1767 + else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360; 5.1768 + } 5.1769 + /* store rounded latitude/longitude values for round-trip safety */ 5.1770 + box->lat_min = pgl_round(lat_min); 5.1771 + box->lat_max = pgl_round(lat_max); 5.1772 + box->lon_min = pgl_round(lon_min); 5.1773 + box->lon_max = pgl_round(lon_max); 5.1774 + /* ensure that rounding does not change orientation */ 5.1775 + if (lon_min > lon_max && box->lon_min == box->lon_max) { 5.1776 + box->lon_min = -180; 5.1777 + box->lon_max = 180; 5.1778 + } 5.1779 +} 5.1780 + 5.1781 +/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */ 5.1782 +PG_FUNCTION_INFO_V1(pgl_create_ebox); 5.1783 +Datum pgl_create_ebox(PG_FUNCTION_ARGS) { 5.1784 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 5.1785 + pgl_ebox_set_boundaries( 5.1786 + box, 5.1787 + PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), 5.1788 + PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3) 5.1789 + ); 5.1790 + PG_RETURN_POINTER(box); 5.1791 +} 5.1792 + 5.1793 +/* create box ("ebox" in SQL) from two points ("epoint"s) */ 5.1794 +/* (can not be used to cover a longitude range of more than 120 degrees) */ 5.1795 +PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints); 5.1796 +Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) { 5.1797 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 5.1798 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 5.1799 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 5.1800 + double lat_min, lat_max, lon_min, lon_max; 5.1801 + double dlon; /* longitude range (delta longitude) */ 5.1802 + /* order latitude and longitude boundaries */ 5.1803 + if (point2->lat < point1->lat) { 5.1804 + lat_min = point2->lat; 5.1805 + lat_max = point1->lat; 5.1806 + } else { 5.1807 + lat_min = point1->lat; 5.1808 + lat_max = point2->lat; 5.1809 + } 5.1810 + if (point2->lon < point1->lon) { 5.1811 + lon_min = point2->lon; 5.1812 + lon_max = point1->lon; 5.1813 + } else { 5.1814 + lon_min = point1->lon; 5.1815 + lon_max = point2->lon; 5.1816 + } 5.1817 + /* calculate longitude range (round to avoid floating point errors) */ 5.1818 + dlon = pgl_round(lon_max - lon_min); 5.1819 + /* determine east-west direction */ 5.1820 + if (dlon >= 240) { 5.1821 + /* assume that 180th meridian is crossed and swap min/max longitude */ 5.1822 + double swap = lon_min; lon_min = lon_max; lon_max = swap; 5.1823 + } else if (dlon > 120) { 5.1824 + /* unclear orientation since delta longitude > 120 */ 5.1825 + ereport(ERROR, ( 5.1826 + errcode(ERRCODE_DATA_EXCEPTION), 5.1827 + errmsg("can not determine east/west orientation for ebox") 5.1828 + )); 5.1829 + } 5.1830 + /* use boundaries to setup box (and perform checks) */ 5.1831 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 5.1832 + /* return result */ 5.1833 + PG_RETURN_POINTER(box); 5.1834 +} 5.1835 + 5.1836 +/* parse box ("ebox" in SQL) */ 5.1837 +/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>' 5.1838 + or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */ 5.1839 +PG_FUNCTION_INFO_V1(pgl_ebox_in); 5.1840 +Datum pgl_ebox_in(PG_FUNCTION_ARGS) { 5.1841 + char *str = PG_GETARG_CSTRING(0); /* input string */ 5.1842 + char *str_lower; /* lower case version of input string */ 5.1843 + char *strptr; /* current position within string */ 5.1844 + int valid; /* number of valid chars */ 5.1845 + int done; /* specifies if latitude or longitude was read */ 5.1846 + double val; /* temporary variable */ 5.1847 + int lat_count = 0; /* count of latitude values parsed */ 5.1848 + int lon_count = 0; /* count of longitufde values parsed */ 5.1849 + double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */ 5.1850 + pgl_box *box; /* return value (to be palloc'ed) */ 5.1851 + /* lowercase input */ 5.1852 + str_lower = psprintf("%s", str); 5.1853 + for (strptr=str_lower; *strptr; strptr++) { 5.1854 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 5.1855 + } 5.1856 + /* reset reading position to start of (lowercase) string */ 5.1857 + strptr = str_lower; 5.1858 + /* check if empty box */ 5.1859 + valid = 0; 5.1860 + sscanf(strptr, " empty %n", &valid); 5.1861 + if (valid && strptr[valid] == 0) { 5.1862 + /* allocate and return empty box */ 5.1863 + box = (pgl_box *)palloc(sizeof(pgl_box)); 5.1864 + pgl_box_set_empty(box); 5.1865 + PG_RETURN_POINTER(box); 5.1866 + } 5.1867 + /* demand four blocks separated by whitespace */ 5.1868 + valid = 0; 5.1869 + sscanf(strptr, " %*s %*s %*s %*s %n", &valid); 5.1870 + /* if four blocks separated by whitespace exist, parse those blocks */ 5.1871 + if (strptr[valid] == 0) while (strptr[0]) { 5.1872 + /* parse either latitude or longitude (whichever found in input string) */ 5.1873 + done = pgl_scan(&strptr, &val, &val); 5.1874 + /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */ 5.1875 + if (done == PGL_SCAN_LAT) { 5.1876 + if (!lat_count) lat_min = val; else lat_max = val; 5.1877 + lat_count++; 5.1878 + } else if (done == PGL_SCAN_LON) { 5.1879 + if (!lon_count) lon_min = val; else lon_max = val; 5.1880 + lon_count++; 5.1881 + } else { 5.1882 + break; 5.1883 + } 5.1884 + } 5.1885 + /* require end of string, and two latitude and two longitude values */ 5.1886 + if (strptr[0] || lat_count != 2 || lon_count != 2) { 5.1887 + ereport(ERROR, ( 5.1888 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.1889 + errmsg("invalid input syntax for type ebox: \"%s\"", str) 5.1890 + )); 5.1891 + } 5.1892 + /* free lower case string */ 5.1893 + pfree(str_lower); 5.1894 + /* order boundaries (maximum greater than minimum) */ 5.1895 + if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; } 5.1896 + if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; } 5.1897 + /* allocate memory for result */ 5.1898 + box = (pgl_box *)palloc(sizeof(pgl_box)); 5.1899 + /* set boundaries (and perform checks) */ 5.1900 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 5.1901 + /* return result */ 5.1902 + PG_RETURN_POINTER(box); 5.1903 +} 5.1904 + 5.1905 +/* set circle to given latitude, longitude, and radius (including checks) */ 5.1906 +static void pgl_ecircle_set_latlon_radius( 5.1907 + pgl_circle *circle, double lat, double lon, double radius 5.1908 +) { 5.1909 + /* set center point (including checks) */ 5.1910 + pgl_epoint_set_latlon(&(circle->center), lat, lon); 5.1911 + /* handle non-positive radius */ 5.1912 + if (isnan(radius)) { 5.1913 + ereport(ERROR, ( 5.1914 + errcode(ERRCODE_DATA_EXCEPTION), 5.1915 + errmsg("invalid radius for ecircle") 5.1916 + )); 5.1917 + } 5.1918 + if (radius == 0) radius = 0; /* avoids -0 */ 5.1919 + else if (radius < 0) { 5.1920 + if (isfinite(radius)) { 5.1921 + ereport(NOTICE, (errmsg("negative radius converted to minus infinity"))); 5.1922 + } 5.1923 + radius = -INFINITY; 5.1924 + } 5.1925 + /* store radius (round-trip safety is ensured by pgl_print_float) */ 5.1926 + circle->radius = radius; 5.1927 +} 5.1928 + 5.1929 +/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */ 5.1930 +PG_FUNCTION_INFO_V1(pgl_create_ecircle); 5.1931 +Datum pgl_create_ecircle(PG_FUNCTION_ARGS) { 5.1932 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 5.1933 + pgl_ecircle_set_latlon_radius( 5.1934 + circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2) 5.1935 + ); 5.1936 + PG_RETURN_POINTER(circle); 5.1937 +} 5.1938 + 5.1939 +/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */ 5.1940 +PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint); 5.1941 +Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) { 5.1942 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.1943 + double radius = PG_GETARG_FLOAT8(1); 5.1944 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 5.1945 + /* set latitude, longitude, radius (and perform checks) */ 5.1946 + pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius); 5.1947 + /* return result */ 5.1948 + PG_RETURN_POINTER(circle); 5.1949 +} 5.1950 + 5.1951 +/* parse circle ("ecircle" in SQL) */ 5.1952 +/* format: '[NS]<float> [EW]<float> <float>' */ 5.1953 +PG_FUNCTION_INFO_V1(pgl_ecircle_in); 5.1954 +Datum pgl_ecircle_in(PG_FUNCTION_ARGS) { 5.1955 + char *str = PG_GETARG_CSTRING(0); /* input string */ 5.1956 + char *strptr = str; /* current position within string */ 5.1957 + double lat, lon, radius; /* parsed values as double precision flaots */ 5.1958 + int valid = 0; /* number of valid chars */ 5.1959 + int done = 0; /* stores if latitude and/or longitude was read */ 5.1960 + pgl_circle *circle; /* return value (to be palloc'ed) */ 5.1961 + /* demand three blocks separated by whitespace */ 5.1962 + sscanf(strptr, " %*s %*s %*s %n", &valid); 5.1963 + /* if three blocks separated by whitespace exist, parse those blocks */ 5.1964 + if (strptr[valid] == 0) { 5.1965 + /* parse latitude and longitude */ 5.1966 + done |= pgl_scan(&strptr, &lat, &lon); 5.1967 + done |= pgl_scan(&strptr, &lat, &lon); 5.1968 + /* parse radius (while incrementing strptr by number of bytes parsed) */ 5.1969 + valid = 0; 5.1970 + if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid; 5.1971 + } 5.1972 + /* require end of string and both latitude and longitude being parsed */ 5.1973 + if (strptr[0] || done != PGL_SCAN_LATLON) { 5.1974 + ereport(ERROR, ( 5.1975 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.1976 + errmsg("invalid input syntax for type ecircle: \"%s\"", str) 5.1977 + )); 5.1978 + } 5.1979 + /* allocate memory for result */ 5.1980 + circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 5.1981 + /* set latitude, longitude, radius (and perform checks) */ 5.1982 + pgl_ecircle_set_latlon_radius(circle, lat, lon, radius); 5.1983 + /* return result */ 5.1984 + PG_RETURN_POINTER(circle); 5.1985 +} 5.1986 + 5.1987 +/* parse cluster ("ecluster" in SQL) */ 5.1988 +PG_FUNCTION_INFO_V1(pgl_ecluster_in); 5.1989 +Datum pgl_ecluster_in(PG_FUNCTION_ARGS) { 5.1990 + int i; 5.1991 + char *str = PG_GETARG_CSTRING(0); /* input string */ 5.1992 + char *str_lower; /* lower case version of input string */ 5.1993 + char *strptr; /* pointer to current reading position of input */ 5.1994 + int npoints_total = 0; /* total number of points in cluster */ 5.1995 + int nentries = 0; /* total number of entries */ 5.1996 + pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */ 5.1997 + int entries_buflen = 4; /* maximum number of elements in entries array */ 5.1998 + int valid; /* number of valid chars processed */ 5.1999 + double lat, lon; /* latitude and longitude of parsed point */ 5.2000 + int entrytype; /* current entry type */ 5.2001 + int npoints; /* number of points in current entry */ 5.2002 + pgl_point *points; /* array of pgl_point for pgl_newentry */ 5.2003 + int points_buflen; /* maximum number of elements in points array */ 5.2004 + int done; /* return value of pgl_scan function */ 5.2005 + pgl_cluster *cluster; /* created cluster */ 5.2006 + /* lowercase input */ 5.2007 + str_lower = psprintf("%s", str); 5.2008 + for (strptr=str_lower; *strptr; strptr++) { 5.2009 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 5.2010 + } 5.2011 + /* reset reading position to start of (lowercase) string */ 5.2012 + strptr = str_lower; 5.2013 + /* allocate initial buffer for entries */ 5.2014 + entries = palloc(entries_buflen * sizeof(pgl_newentry)); 5.2015 + /* parse until end of string */ 5.2016 + while (strptr[0]) { 5.2017 + /* require previous white-space or closing parenthesis before next token */ 5.2018 + if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') { 5.2019 + goto pgl_ecluster_in_error; 5.2020 + } 5.2021 + /* ignore token "empty" */ 5.2022 + valid = 0; sscanf(strptr, " empty %n", &valid); 5.2023 + if (valid) { strptr += valid; continue; } 5.2024 + /* test for "point" token */ 5.2025 + valid = 0; sscanf(strptr, " point ( %n", &valid); 5.2026 + if (valid) { 5.2027 + strptr += valid; 5.2028 + entrytype = PGL_ENTRY_POINT; 5.2029 + goto pgl_ecluster_in_type_ok; 5.2030 + } 5.2031 + /* test for "path" token */ 5.2032 + valid = 0; sscanf(strptr, " path ( %n", &valid); 5.2033 + if (valid) { 5.2034 + strptr += valid; 5.2035 + entrytype = PGL_ENTRY_PATH; 5.2036 + goto pgl_ecluster_in_type_ok; 5.2037 + } 5.2038 + /* test for "outline" token */ 5.2039 + valid = 0; sscanf(strptr, " outline ( %n", &valid); 5.2040 + if (valid) { 5.2041 + strptr += valid; 5.2042 + entrytype = PGL_ENTRY_OUTLINE; 5.2043 + goto pgl_ecluster_in_type_ok; 5.2044 + } 5.2045 + /* test for "polygon" token */ 5.2046 + valid = 0; sscanf(strptr, " polygon ( %n", &valid); 5.2047 + if (valid) { 5.2048 + strptr += valid; 5.2049 + entrytype = PGL_ENTRY_POLYGON; 5.2050 + goto pgl_ecluster_in_type_ok; 5.2051 + } 5.2052 + /* error if no valid token found */ 5.2053 + goto pgl_ecluster_in_error; 5.2054 + pgl_ecluster_in_type_ok: 5.2055 + /* check if pgl_newentry array needs to grow */ 5.2056 + if (nentries == entries_buflen) { 5.2057 + pgl_newentry *newbuf; 5.2058 + entries_buflen *= 2; 5.2059 + newbuf = palloc(entries_buflen * sizeof(pgl_newentry)); 5.2060 + memcpy(newbuf, entries, nentries * sizeof(pgl_newentry)); 5.2061 + pfree(entries); 5.2062 + entries = newbuf; 5.2063 + } 5.2064 + /* reset number of points for current entry */ 5.2065 + npoints = 0; 5.2066 + /* allocate array for points */ 5.2067 + points_buflen = 4; 5.2068 + points = palloc(points_buflen * sizeof(pgl_point)); 5.2069 + /* parse until closing parenthesis */ 5.2070 + while (strptr[0] != ')') { 5.2071 + /* error on unexpected end of string */ 5.2072 + if (strptr[0] == 0) goto pgl_ecluster_in_error; 5.2073 + /* mark neither latitude nor longitude as read */ 5.2074 + done = PGL_SCAN_NONE; 5.2075 + /* require white-space before second, third, etc. point */ 5.2076 + if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 5.2077 + /* scan latitude (or longitude) */ 5.2078 + done |= pgl_scan(&strptr, &lat, &lon); 5.2079 + /* require white-space before second coordinate */ 5.2080 + if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 5.2081 + /* scan longitude (or latitude) */ 5.2082 + done |= pgl_scan(&strptr, &lat, &lon); 5.2083 + /* error unless both latitude and longitude were parsed */ 5.2084 + if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error; 5.2085 + /* throw error if number of points is too high */ 5.2086 + if (npoints_total == PGL_CLUSTER_MAXPOINTS) { 5.2087 + ereport(ERROR, ( 5.2088 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.2089 + errmsg( 5.2090 + "too many points for ecluster entry (maximum %i)", 5.2091 + PGL_CLUSTER_MAXPOINTS 5.2092 + ) 5.2093 + )); 5.2094 + } 5.2095 + /* check if pgl_point array needs to grow */ 5.2096 + if (npoints == points_buflen) { 5.2097 + pgl_point *newbuf; 5.2098 + points_buflen *= 2; 5.2099 + newbuf = palloc(points_buflen * sizeof(pgl_point)); 5.2100 + memcpy(newbuf, points, npoints * sizeof(pgl_point)); 5.2101 + pfree(points); 5.2102 + points = newbuf; 5.2103 + } 5.2104 + /* append point to pgl_point array (includes checks) */ 5.2105 + pgl_epoint_set_latlon(&(points[npoints++]), lat, lon); 5.2106 + /* increase total number of points */ 5.2107 + npoints_total++; 5.2108 + } 5.2109 + /* error if entry has no points */ 5.2110 + if (!npoints) goto pgl_ecluster_in_error; 5.2111 + /* entries with one point are automatically of type "point" */ 5.2112 + if (npoints == 1) entrytype = PGL_ENTRY_POINT; 5.2113 + /* if entries have more than one point */ 5.2114 + else { 5.2115 + /* throw error if entry type is "point" */ 5.2116 + if (entrytype == PGL_ENTRY_POINT) { 5.2117 + ereport(ERROR, ( 5.2118 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.2119 + errmsg("invalid input syntax for type ecluster (point entry with more than one point)") 5.2120 + )); 5.2121 + } 5.2122 + /* coerce outlines and polygons with more than 2 points to be a path */ 5.2123 + if (npoints == 2) entrytype = PGL_ENTRY_PATH; 5.2124 + } 5.2125 + /* append entry to pgl_newentry array */ 5.2126 + entries[nentries].entrytype = entrytype; 5.2127 + entries[nentries].npoints = npoints; 5.2128 + entries[nentries].points = points; 5.2129 + nentries++; 5.2130 + /* consume closing parenthesis */ 5.2131 + strptr++; 5.2132 + /* consume white-space */ 5.2133 + while (isspace(strptr[0])) strptr++; 5.2134 + } 5.2135 + /* free lower case string */ 5.2136 + pfree(str_lower); 5.2137 + /* create cluster from pgl_newentry array */ 5.2138 + cluster = pgl_new_cluster(nentries, entries); 5.2139 + /* free pgl_newentry array */ 5.2140 + for (i=0; i<nentries; i++) pfree(entries[i].points); 5.2141 + pfree(entries); 5.2142 + /* set bounding circle of cluster and check east/west orientation */ 5.2143 + if (!pgl_finalize_cluster(cluster)) { 5.2144 + ereport(ERROR, ( 5.2145 + errcode(ERRCODE_DATA_EXCEPTION), 5.2146 + errmsg("can not determine east/west orientation for ecluster"), 5.2147 + errhint("Ensure that each entry has a longitude span of less than 180 degrees.") 5.2148 + )); 5.2149 + } 5.2150 + /* return cluster */ 5.2151 + PG_RETURN_POINTER(cluster); 5.2152 + /* code to throw error */ 5.2153 + pgl_ecluster_in_error: 5.2154 + ereport(ERROR, ( 5.2155 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 5.2156 + errmsg("invalid input syntax for type ecluster: \"%s\"", str) 5.2157 + )); 5.2158 +} 5.2159 + 5.2160 +/* convert point ("epoint") to string representation */ 5.2161 +PG_FUNCTION_INFO_V1(pgl_epoint_out); 5.2162 +Datum pgl_epoint_out(PG_FUNCTION_ARGS) { 5.2163 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2164 + char latstr[PGL_NUMBUFLEN]; 5.2165 + char lonstr[PGL_NUMBUFLEN]; 5.2166 + pgl_print_lat(latstr, point->lat); 5.2167 + pgl_print_lon(lonstr, point->lon); 5.2168 + PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr)); 5.2169 +} 5.2170 + 5.2171 +/* convert box ("ebox") to string representation */ 5.2172 +PG_FUNCTION_INFO_V1(pgl_ebox_out); 5.2173 +Datum pgl_ebox_out(PG_FUNCTION_ARGS) { 5.2174 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 5.2175 + double lon_min = box->lon_min; 5.2176 + double lon_max = box->lon_max; 5.2177 + char lat_min_str[PGL_NUMBUFLEN]; 5.2178 + char lat_max_str[PGL_NUMBUFLEN]; 5.2179 + char lon_min_str[PGL_NUMBUFLEN]; 5.2180 + char lon_max_str[PGL_NUMBUFLEN]; 5.2181 + /* return string "empty" if box is set to be empty */ 5.2182 + if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty"); 5.2183 + /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */ 5.2184 + /* (required since pgl_box_in orders the longitude boundaries) */ 5.2185 + if (lon_min > lon_max) { 5.2186 + if (lon_min + lon_max >= 0) lon_min -= 360; 5.2187 + else lon_max += 360; 5.2188 + } 5.2189 + /* format and return result */ 5.2190 + pgl_print_lat(lat_min_str, box->lat_min); 5.2191 + pgl_print_lat(lat_max_str, box->lat_max); 5.2192 + pgl_print_lon(lon_min_str, lon_min); 5.2193 + pgl_print_lon(lon_max_str, lon_max); 5.2194 + PG_RETURN_CSTRING(psprintf( 5.2195 + "%s %s %s %s", 5.2196 + lat_min_str, lon_min_str, lat_max_str, lon_max_str 5.2197 + )); 5.2198 +} 5.2199 + 5.2200 +/* convert circle ("ecircle") to string representation */ 5.2201 +PG_FUNCTION_INFO_V1(pgl_ecircle_out); 5.2202 +Datum pgl_ecircle_out(PG_FUNCTION_ARGS) { 5.2203 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 5.2204 + char latstr[PGL_NUMBUFLEN]; 5.2205 + char lonstr[PGL_NUMBUFLEN]; 5.2206 + char radstr[PGL_NUMBUFLEN]; 5.2207 + pgl_print_lat(latstr, circle->center.lat); 5.2208 + pgl_print_lon(lonstr, circle->center.lon); 5.2209 + pgl_print_float(radstr, circle->radius); 5.2210 + PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr)); 5.2211 +} 5.2212 + 5.2213 +/* convert cluster ("ecluster") to string representation */ 5.2214 +PG_FUNCTION_INFO_V1(pgl_ecluster_out); 5.2215 +Datum pgl_ecluster_out(PG_FUNCTION_ARGS) { 5.2216 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2217 + char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */ 5.2218 + char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */ 5.2219 + char ***strings; /* array of array of strings */ 5.2220 + char *string; /* string of current token */ 5.2221 + char *res, *resptr; /* result and pointer to current write position */ 5.2222 + size_t reslen = 1; /* length of result (init with 1 for terminator) */ 5.2223 + int npoints; /* number of points of current entry */ 5.2224 + int i, j; /* i: entry, j: point in entry */ 5.2225 + /* handle empty clusters */ 5.2226 + if (cluster->nentries == 0) { 5.2227 + /* free detoasted cluster (if copy) */ 5.2228 + PG_FREE_IF_COPY(cluster, 0); 5.2229 + /* return static result */ 5.2230 + PG_RETURN_CSTRING("empty"); 5.2231 + } 5.2232 + /* allocate array of array of strings */ 5.2233 + strings = palloc(cluster->nentries * sizeof(char **)); 5.2234 + /* iterate over all entries in cluster */ 5.2235 + for (i=0; i<cluster->nentries; i++) { 5.2236 + /* get number of points in entry */ 5.2237 + npoints = cluster->entries[i].npoints; 5.2238 + /* allocate array of strings (one string for each point plus two extra) */ 5.2239 + strings[i] = palloc((2 + npoints) * sizeof(char *)); 5.2240 + /* determine opening string */ 5.2241 + switch (cluster->entries[i].entrytype) { 5.2242 + case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break; 5.2243 + case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break; 5.2244 + case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break; 5.2245 + case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break; 5.2246 + default: string = (i==0)?"unknown" :" unknown"; 5.2247 + } 5.2248 + /* use opening string as first string in array */ 5.2249 + strings[i][0] = string; 5.2250 + /* update result length (for allocating result string later) */ 5.2251 + reslen += strlen(string); 5.2252 + /* iterate over all points */ 5.2253 + for (j=0; j<npoints; j++) { 5.2254 + /* create string representation of point */ 5.2255 + pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat); 5.2256 + pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon); 5.2257 + string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr); 5.2258 + /* copy string pointer to string array */ 5.2259 + strings[i][j+1] = string; 5.2260 + /* update result length (for allocating result string later) */ 5.2261 + reslen += strlen(string); 5.2262 + } 5.2263 + /* use closing parenthesis as last string in array */ 5.2264 + strings[i][npoints+1] = ")"; 5.2265 + /* update result length (for allocating result string later) */ 5.2266 + reslen++; 5.2267 + } 5.2268 + /* allocate result string */ 5.2269 + res = palloc(reslen); 5.2270 + /* set write pointer to begin of result string */ 5.2271 + resptr = res; 5.2272 + /* copy strings into result string */ 5.2273 + for (i=0; i<cluster->nentries; i++) { 5.2274 + npoints = cluster->entries[i].npoints; 5.2275 + for (j=0; j<npoints+2; j++) { 5.2276 + string = strings[i][j]; 5.2277 + strcpy(resptr, string); 5.2278 + resptr += strlen(string); 5.2279 + /* free strings allocated by psprintf */ 5.2280 + if (j != 0 && j != npoints+1) pfree(string); 5.2281 + } 5.2282 + /* free array of strings */ 5.2283 + pfree(strings[i]); 5.2284 + } 5.2285 + /* free array of array of strings */ 5.2286 + pfree(strings); 5.2287 + /* free detoasted cluster (if copy) */ 5.2288 + PG_FREE_IF_COPY(cluster, 0); 5.2289 + /* return result */ 5.2290 + PG_RETURN_CSTRING(res); 5.2291 +} 5.2292 + 5.2293 +/* binary input function for point ("epoint") */ 5.2294 +PG_FUNCTION_INFO_V1(pgl_epoint_recv); 5.2295 +Datum pgl_epoint_recv(PG_FUNCTION_ARGS) { 5.2296 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 5.2297 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 5.2298 + point->lat = pq_getmsgfloat8(buf); 5.2299 + point->lon = pq_getmsgfloat8(buf); 5.2300 + PG_RETURN_POINTER(point); 5.2301 +} 5.2302 + 5.2303 +/* binary input function for box ("ebox") */ 5.2304 +PG_FUNCTION_INFO_V1(pgl_ebox_recv); 5.2305 +Datum pgl_ebox_recv(PG_FUNCTION_ARGS) { 5.2306 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 5.2307 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 5.2308 + box->lat_min = pq_getmsgfloat8(buf); 5.2309 + box->lat_max = pq_getmsgfloat8(buf); 5.2310 + box->lon_min = pq_getmsgfloat8(buf); 5.2311 + box->lon_max = pq_getmsgfloat8(buf); 5.2312 + PG_RETURN_POINTER(box); 5.2313 +} 5.2314 + 5.2315 +/* binary input function for circle ("ecircle") */ 5.2316 +PG_FUNCTION_INFO_V1(pgl_ecircle_recv); 5.2317 +Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) { 5.2318 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 5.2319 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 5.2320 + circle->center.lat = pq_getmsgfloat8(buf); 5.2321 + circle->center.lon = pq_getmsgfloat8(buf); 5.2322 + circle->radius = pq_getmsgfloat8(buf); 5.2323 + PG_RETURN_POINTER(circle); 5.2324 +} 5.2325 + 5.2326 +/* TODO: binary receive function for cluster */ 5.2327 + 5.2328 +/* binary output function for point ("epoint") */ 5.2329 +PG_FUNCTION_INFO_V1(pgl_epoint_send); 5.2330 +Datum pgl_epoint_send(PG_FUNCTION_ARGS) { 5.2331 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2332 + StringInfoData buf; 5.2333 + pq_begintypsend(&buf); 5.2334 + pq_sendfloat8(&buf, point->lat); 5.2335 + pq_sendfloat8(&buf, point->lon); 5.2336 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 5.2337 +} 5.2338 + 5.2339 +/* binary output function for box ("ebox") */ 5.2340 +PG_FUNCTION_INFO_V1(pgl_ebox_send); 5.2341 +Datum pgl_ebox_send(PG_FUNCTION_ARGS) { 5.2342 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 5.2343 + StringInfoData buf; 5.2344 + pq_begintypsend(&buf); 5.2345 + pq_sendfloat8(&buf, box->lat_min); 5.2346 + pq_sendfloat8(&buf, box->lat_max); 5.2347 + pq_sendfloat8(&buf, box->lon_min); 5.2348 + pq_sendfloat8(&buf, box->lon_max); 5.2349 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 5.2350 +} 5.2351 + 5.2352 +/* binary output function for circle ("ecircle") */ 5.2353 +PG_FUNCTION_INFO_V1(pgl_ecircle_send); 5.2354 +Datum pgl_ecircle_send(PG_FUNCTION_ARGS) { 5.2355 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 5.2356 + StringInfoData buf; 5.2357 + pq_begintypsend(&buf); 5.2358 + pq_sendfloat8(&buf, circle->center.lat); 5.2359 + pq_sendfloat8(&buf, circle->center.lon); 5.2360 + pq_sendfloat8(&buf, circle->radius); 5.2361 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 5.2362 +} 5.2363 + 5.2364 +/* TODO: binary send functions for cluster */ 5.2365 + 5.2366 +/* cast point ("epoint") to box ("ebox") */ 5.2367 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox); 5.2368 +Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) { 5.2369 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2370 + pgl_box *box = palloc(sizeof(pgl_box)); 5.2371 + box->lat_min = point->lat; 5.2372 + box->lat_max = point->lat; 5.2373 + box->lon_min = point->lon; 5.2374 + box->lon_max = point->lon; 5.2375 + PG_RETURN_POINTER(box); 5.2376 +} 5.2377 + 5.2378 +/* cast point ("epoint") to circle ("ecircle") */ 5.2379 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle); 5.2380 +Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) { 5.2381 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2382 + pgl_circle *circle = palloc(sizeof(pgl_box)); 5.2383 + circle->center = *point; 5.2384 + circle->radius = 0; 5.2385 + PG_RETURN_POINTER(circle); 5.2386 +} 5.2387 + 5.2388 +/* cast point ("epoint") to cluster ("ecluster") */ 5.2389 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster); 5.2390 +Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) { 5.2391 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2392 + pgl_newentry entry; 5.2393 + pgl_cluster *cluster; 5.2394 + entry.entrytype = PGL_ENTRY_POINT; 5.2395 + entry.npoints = 1; 5.2396 + entry.points = point; 5.2397 + cluster = pgl_new_cluster(1, &entry); 5.2398 + pgl_finalize_cluster(cluster); /* NOTE: should not fail */ 5.2399 + PG_RETURN_POINTER(cluster); 5.2400 +} 5.2401 + 5.2402 +/* cast box ("ebox") to cluster ("ecluster") */ 5.2403 +#define pgl_ebox_to_ecluster_macro(i, a, b) \ 5.2404 + entries[i].entrytype = PGL_ENTRY_POLYGON; \ 5.2405 + entries[i].npoints = 4; \ 5.2406 + entries[i].points = points[i]; \ 5.2407 + points[i][0].lat = box->lat_min; \ 5.2408 + points[i][0].lon = (a); \ 5.2409 + points[i][1].lat = box->lat_min; \ 5.2410 + points[i][1].lon = (b); \ 5.2411 + points[i][2].lat = box->lat_max; \ 5.2412 + points[i][2].lon = (b); \ 5.2413 + points[i][3].lat = box->lat_max; \ 5.2414 + points[i][3].lon = (a); 5.2415 +PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster); 5.2416 +Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) { 5.2417 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 5.2418 + double lon, dlon; 5.2419 + int nentries; 5.2420 + pgl_newentry entries[3]; 5.2421 + pgl_point points[3][4]; 5.2422 + pgl_cluster *cluster; 5.2423 + if (box->lat_min > box->lat_max) { 5.2424 + nentries = 0; 5.2425 + } else if (box->lon_min > box->lon_max) { 5.2426 + if (box->lon_min < 0) { 5.2427 + lon = pgl_round((box->lon_min + 180) / 2.0); 5.2428 + nentries = 3; 5.2429 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 5.2430 + pgl_ebox_to_ecluster_macro(1, lon, 180); 5.2431 + pgl_ebox_to_ecluster_macro(2, -180, box->lon_max); 5.2432 + } else if (box->lon_max > 0) { 5.2433 + lon = pgl_round((box->lon_max - 180) / 2.0); 5.2434 + nentries = 3; 5.2435 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 5.2436 + pgl_ebox_to_ecluster_macro(1, -180, lon); 5.2437 + pgl_ebox_to_ecluster_macro(2, lon, box->lon_max); 5.2438 + } else { 5.2439 + nentries = 2; 5.2440 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 5.2441 + pgl_ebox_to_ecluster_macro(1, -180, box->lon_max); 5.2442 + } 5.2443 + } else { 5.2444 + dlon = pgl_round(box->lon_max - box->lon_min); 5.2445 + if (dlon < 180) { 5.2446 + nentries = 1; 5.2447 + pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max); 5.2448 + } else { 5.2449 + lon = pgl_round((box->lon_min + box->lon_max) / 2.0); 5.2450 + if ( 5.2451 + pgl_round(lon - box->lon_min) < 180 && 5.2452 + pgl_round(box->lon_max - lon) < 180 5.2453 + ) { 5.2454 + nentries = 2; 5.2455 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 5.2456 + pgl_ebox_to_ecluster_macro(1, lon, box->lon_max); 5.2457 + } else { 5.2458 + nentries = 3; 5.2459 + pgl_ebox_to_ecluster_macro(0, box->lon_min, -60); 5.2460 + pgl_ebox_to_ecluster_macro(1, -60, 60); 5.2461 + pgl_ebox_to_ecluster_macro(2, 60, box->lon_max); 5.2462 + } 5.2463 + } 5.2464 + } 5.2465 + cluster = pgl_new_cluster(nentries, entries); 5.2466 + pgl_finalize_cluster(cluster); /* NOTE: should not fail */ 5.2467 + PG_RETURN_POINTER(cluster); 5.2468 +} 5.2469 + 5.2470 +/* extract latitude from point ("epoint") */ 5.2471 +PG_FUNCTION_INFO_V1(pgl_epoint_lat); 5.2472 +Datum pgl_epoint_lat(PG_FUNCTION_ARGS) { 5.2473 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat); 5.2474 +} 5.2475 + 5.2476 +/* extract longitude from point ("epoint") */ 5.2477 +PG_FUNCTION_INFO_V1(pgl_epoint_lon); 5.2478 +Datum pgl_epoint_lon(PG_FUNCTION_ARGS) { 5.2479 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon); 5.2480 +} 5.2481 + 5.2482 +/* extract minimum latitude from box ("ebox") */ 5.2483 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_min); 5.2484 +Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) { 5.2485 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min); 5.2486 +} 5.2487 + 5.2488 +/* extract maximum latitude from box ("ebox") */ 5.2489 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_max); 5.2490 +Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) { 5.2491 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max); 5.2492 +} 5.2493 + 5.2494 +/* extract minimum longitude from box ("ebox") */ 5.2495 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_min); 5.2496 +Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) { 5.2497 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min); 5.2498 +} 5.2499 + 5.2500 +/* extract maximum longitude from box ("ebox") */ 5.2501 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_max); 5.2502 +Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) { 5.2503 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max); 5.2504 +} 5.2505 + 5.2506 +/* extract center point from circle ("ecircle") */ 5.2507 +PG_FUNCTION_INFO_V1(pgl_ecircle_center); 5.2508 +Datum pgl_ecircle_center(PG_FUNCTION_ARGS) { 5.2509 + PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center)); 5.2510 +} 5.2511 + 5.2512 +/* extract radius from circle ("ecircle") */ 5.2513 +PG_FUNCTION_INFO_V1(pgl_ecircle_radius); 5.2514 +Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) { 5.2515 + PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius); 5.2516 +} 5.2517 + 5.2518 +/* check if point is inside box (overlap operator "&&") in SQL */ 5.2519 +PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap); 5.2520 +Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) { 5.2521 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2522 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1); 5.2523 + PG_RETURN_BOOL(pgl_point_in_box(point, box)); 5.2524 +} 5.2525 + 5.2526 +/* check if point is inside circle (overlap operator "&&") in SQL */ 5.2527 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap); 5.2528 +Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) { 5.2529 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2530 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 5.2531 + PG_RETURN_BOOL( 5.2532 + pgl_distance( 5.2533 + point->lat, point->lon, 5.2534 + circle->center.lat, circle->center.lon 5.2535 + ) <= circle->radius 5.2536 + ); 5.2537 +} 5.2538 + 5.2539 +/* check if point is inside cluster (overlap operator "&&") in SQL */ 5.2540 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap); 5.2541 +Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) { 5.2542 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2543 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2544 + bool retval; 5.2545 + /* points outside bounding circle are always assumed to be non-overlapping 5.2546 + (necessary for consistent table and index scans) */ 5.2547 + if ( 5.2548 + pgl_distance( 5.2549 + point->lat, point->lon, 5.2550 + cluster->bounding.center.lat, cluster->bounding.center.lon 5.2551 + ) > cluster->bounding.radius 5.2552 + ) retval = false; 5.2553 + else retval = pgl_point_in_cluster(point, cluster, false); 5.2554 + PG_FREE_IF_COPY(cluster, 1); 5.2555 + PG_RETURN_BOOL(retval); 5.2556 +} 5.2557 + 5.2558 +/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */ 5.2559 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap); 5.2560 +Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) { 5.2561 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2562 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2563 + bool retval = pgl_distance( 5.2564 + point->lat, point->lon, 5.2565 + cluster->bounding.center.lat, cluster->bounding.center.lon 5.2566 + ) <= cluster->bounding.radius; 5.2567 + PG_FREE_IF_COPY(cluster, 1); 5.2568 + PG_RETURN_BOOL(retval); 5.2569 +} 5.2570 + 5.2571 +/* check if two boxes overlap (overlap operator "&&") in SQL */ 5.2572 +PG_FUNCTION_INFO_V1(pgl_ebox_overlap); 5.2573 +Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) { 5.2574 + pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0); 5.2575 + pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1); 5.2576 + PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2)); 5.2577 +} 5.2578 + 5.2579 +/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */ 5.2580 +PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap); 5.2581 +Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) { 5.2582 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 5.2583 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 5.2584 + PG_RETURN_BOOL( 5.2585 + pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius 5.2586 + ); 5.2587 +} 5.2588 + 5.2589 +/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */ 5.2590 +PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap); 5.2591 +Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) { 5.2592 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 5.2593 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2594 + bool retval = pgl_estimate_point_box_distance( 5.2595 + &cluster->bounding.center, 5.2596 + box 5.2597 + ) <= cluster->bounding.radius; 5.2598 + PG_FREE_IF_COPY(cluster, 1); 5.2599 + PG_RETURN_BOOL(retval); 5.2600 +} 5.2601 + 5.2602 +/* check if two circles overlap (overlap operator "&&") in SQL */ 5.2603 +PG_FUNCTION_INFO_V1(pgl_ecircle_overlap); 5.2604 +Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) { 5.2605 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 5.2606 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 5.2607 + PG_RETURN_BOOL( 5.2608 + pgl_distance( 5.2609 + circle1->center.lat, circle1->center.lon, 5.2610 + circle2->center.lat, circle2->center.lon 5.2611 + ) <= circle1->radius + circle2->radius 5.2612 + ); 5.2613 +} 5.2614 + 5.2615 +/* check if circle and cluster overlap (overlap operator "&&") in SQL */ 5.2616 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap); 5.2617 +Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) { 5.2618 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 5.2619 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2620 + bool retval = ( 5.2621 + pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius 5.2622 + ); 5.2623 + PG_FREE_IF_COPY(cluster, 1); 5.2624 + PG_RETURN_BOOL(retval); 5.2625 +} 5.2626 + 5.2627 +/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */ 5.2628 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap); 5.2629 +Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) { 5.2630 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 5.2631 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2632 + bool retval = pgl_distance( 5.2633 + circle->center.lat, circle->center.lon, 5.2634 + cluster->bounding.center.lat, cluster->bounding.center.lon 5.2635 + ) <= circle->radius + cluster->bounding.radius; 5.2636 + PG_FREE_IF_COPY(cluster, 1); 5.2637 + PG_RETURN_BOOL(retval); 5.2638 +} 5.2639 + 5.2640 +/* check if two clusters overlap (overlap operator "&&") in SQL */ 5.2641 +PG_FUNCTION_INFO_V1(pgl_ecluster_overlap); 5.2642 +Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) { 5.2643 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2644 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2645 + bool retval; 5.2646 + /* clusters with non-touching bounding circles are always assumed to be 5.2647 + non-overlapping (improves performance and is necessary for consistent 5.2648 + table and index scans) */ 5.2649 + if ( 5.2650 + pgl_distance( 5.2651 + cluster1->bounding.center.lat, cluster1->bounding.center.lon, 5.2652 + cluster2->bounding.center.lat, cluster2->bounding.center.lon 5.2653 + ) > cluster1->bounding.radius + cluster2->bounding.radius 5.2654 + ) retval = false; 5.2655 + else retval = pgl_clusters_overlap(cluster1, cluster2); 5.2656 + PG_FREE_IF_COPY(cluster1, 0); 5.2657 + PG_FREE_IF_COPY(cluster2, 1); 5.2658 + PG_RETURN_BOOL(retval); 5.2659 +} 5.2660 + 5.2661 +/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */ 5.2662 +PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap); 5.2663 +Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) { 5.2664 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2665 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2666 + bool retval = pgl_distance( 5.2667 + cluster1->bounding.center.lat, cluster1->bounding.center.lon, 5.2668 + cluster2->bounding.center.lat, cluster2->bounding.center.lon 5.2669 + ) <= cluster1->bounding.radius + cluster2->bounding.radius; 5.2670 + PG_FREE_IF_COPY(cluster1, 0); 5.2671 + PG_FREE_IF_COPY(cluster2, 1); 5.2672 + PG_RETURN_BOOL(retval); 5.2673 +} 5.2674 + 5.2675 +/* check if second cluster is in first cluster (cont. operator "@>) in SQL */ 5.2676 +PG_FUNCTION_INFO_V1(pgl_ecluster_contains); 5.2677 +Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) { 5.2678 + pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2679 + pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2680 + bool retval; 5.2681 + /* clusters with non-touching bounding circles are always assumed to be 5.2682 + non-overlapping (improves performance and is necessary for consistent 5.2683 + table and index scans) */ 5.2684 + if ( 5.2685 + pgl_distance( 5.2686 + outer->bounding.center.lat, outer->bounding.center.lon, 5.2687 + inner->bounding.center.lat, inner->bounding.center.lon 5.2688 + ) > outer->bounding.radius + inner->bounding.radius 5.2689 + ) retval = false; 5.2690 + else retval = pgl_cluster_in_cluster(outer, inner); 5.2691 + PG_FREE_IF_COPY(outer, 0); 5.2692 + PG_FREE_IF_COPY(inner, 1); 5.2693 + PG_RETURN_BOOL(retval); 5.2694 +} 5.2695 + 5.2696 +/* calculate distance between two points ("<->" operator) in SQL */ 5.2697 +PG_FUNCTION_INFO_V1(pgl_epoint_distance); 5.2698 +Datum pgl_epoint_distance(PG_FUNCTION_ARGS) { 5.2699 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 5.2700 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 5.2701 + PG_RETURN_FLOAT8(pgl_distance( 5.2702 + point1->lat, point1->lon, point2->lat, point2->lon 5.2703 + )); 5.2704 +} 5.2705 + 5.2706 +/* calculate point to circle distance ("<->" operator) in SQL */ 5.2707 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance); 5.2708 +Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) { 5.2709 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2710 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 5.2711 + double distance = pgl_distance( 5.2712 + point->lat, point->lon, circle->center.lat, circle->center.lon 5.2713 + ) - circle->radius; 5.2714 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 5.2715 +} 5.2716 + 5.2717 +/* calculate point to cluster distance ("<->" operator) in SQL */ 5.2718 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance); 5.2719 +Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) { 5.2720 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 5.2721 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2722 + double distance = pgl_point_cluster_distance(point, cluster); 5.2723 + PG_FREE_IF_COPY(cluster, 1); 5.2724 + PG_RETURN_FLOAT8(distance); 5.2725 +} 5.2726 + 5.2727 +/* calculate distance between two circles ("<->" operator) in SQL */ 5.2728 +PG_FUNCTION_INFO_V1(pgl_ecircle_distance); 5.2729 +Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) { 5.2730 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 5.2731 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 5.2732 + double distance = pgl_distance( 5.2733 + circle1->center.lat, circle1->center.lon, 5.2734 + circle2->center.lat, circle2->center.lon 5.2735 + ) - (circle1->radius + circle2->radius); 5.2736 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 5.2737 +} 5.2738 + 5.2739 +/* calculate circle to cluster distance ("<->" operator) in SQL */ 5.2740 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance); 5.2741 +Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) { 5.2742 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 5.2743 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2744 + double distance = ( 5.2745 + pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius 5.2746 + ); 5.2747 + PG_FREE_IF_COPY(cluster, 1); 5.2748 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 5.2749 +} 5.2750 + 5.2751 +/* calculate distance between two clusters ("<->" operator) in SQL */ 5.2752 +PG_FUNCTION_INFO_V1(pgl_ecluster_distance); 5.2753 +Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) { 5.2754 + pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2755 + pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2756 + double retval = pgl_cluster_distance(cluster1, cluster2); 5.2757 + PG_FREE_IF_COPY(cluster1, 0); 5.2758 + PG_FREE_IF_COPY(cluster2, 1); 5.2759 + PG_RETURN_FLOAT8(retval); 5.2760 +} 5.2761 + 5.2762 +PG_FUNCTION_INFO_V1(pgl_ecluster_monte_carlo_area); 5.2763 +Datum pgl_ecluster_monte_carlo_area(PG_FUNCTION_ARGS) { 5.2764 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2765 + int samples = PG_GETARG_INT32(1); 5.2766 + double retval = pgl_monte_carlo_area(cluster, samples); 5.2767 + PG_FREE_IF_COPY(cluster, 0); 5.2768 + PG_RETURN_FLOAT8(retval); 5.2769 +} 5.2770 + 5.2771 +PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_fair_distance); 5.2772 +Datum pgl_ecluster_epoint_fair_distance(PG_FUNCTION_ARGS) { 5.2773 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 5.2774 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(1); 5.2775 + int samples = PG_GETARG_INT32(2); 5.2776 + double retval = pgl_fair_distance(point, cluster, samples); 5.2777 + PG_FREE_IF_COPY(cluster, 0); 5.2778 + PG_RETURN_FLOAT8(retval); 5.2779 +} 5.2780 + 5.2781 + 5.2782 +/*-----------------------------------------------------------* 5.2783 + * B-tree comparison operators and index support functions * 5.2784 + *-----------------------------------------------------------*/ 5.2785 + 5.2786 +/* macro for a B-tree operator (without detoasting) */ 5.2787 +#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \ 5.2788 + PG_FUNCTION_INFO_V1(func); \ 5.2789 + Datum func(PG_FUNCTION_ARGS) { \ 5.2790 + type *a = (type *)PG_GETARG_POINTER(0); \ 5.2791 + type *b = (type *)PG_GETARG_POINTER(1); \ 5.2792 + PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \ 5.2793 + } 5.2794 + 5.2795 +/* macro for a B-tree comparison function (without detoasting) */ 5.2796 +#define PGL_BTREE_CMP(func, type, cmpfunc) \ 5.2797 + PG_FUNCTION_INFO_V1(func); \ 5.2798 + Datum func(PG_FUNCTION_ARGS) { \ 5.2799 + type *a = (type *)PG_GETARG_POINTER(0); \ 5.2800 + type *b = (type *)PG_GETARG_POINTER(1); \ 5.2801 + PG_RETURN_INT32(cmpfunc(a, b)); \ 5.2802 + } 5.2803 + 5.2804 +/* macro for a B-tree operator (with detoasting) */ 5.2805 +#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \ 5.2806 + PG_FUNCTION_INFO_V1(func); \ 5.2807 + Datum func(PG_FUNCTION_ARGS) { \ 5.2808 + bool res; \ 5.2809 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 5.2810 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 5.2811 + res = cmpfunc(a, b) oper 0; \ 5.2812 + PG_FREE_IF_COPY(a, 0); \ 5.2813 + PG_FREE_IF_COPY(b, 1); \ 5.2814 + PG_RETURN_BOOL(res); \ 5.2815 + } 5.2816 + 5.2817 +/* macro for a B-tree comparison function (with detoasting) */ 5.2818 +#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \ 5.2819 + PG_FUNCTION_INFO_V1(func); \ 5.2820 + Datum func(PG_FUNCTION_ARGS) { \ 5.2821 + int32_t res; \ 5.2822 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 5.2823 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 5.2824 + res = cmpfunc(a, b); \ 5.2825 + PG_FREE_IF_COPY(a, 0); \ 5.2826 + PG_FREE_IF_COPY(b, 1); \ 5.2827 + PG_RETURN_INT32(res); \ 5.2828 + } 5.2829 + 5.2830 +/* B-tree operators and comparison function for point */ 5.2831 +PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <) 5.2832 +PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=) 5.2833 +PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==) 5.2834 +PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=) 5.2835 +PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=) 5.2836 +PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >) 5.2837 +PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp) 5.2838 + 5.2839 +/* B-tree operators and comparison function for box */ 5.2840 +PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <) 5.2841 +PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=) 5.2842 +PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==) 5.2843 +PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=) 5.2844 +PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=) 5.2845 +PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >) 5.2846 +PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp) 5.2847 + 5.2848 +/* B-tree operators and comparison function for circle */ 5.2849 +PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <) 5.2850 +PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=) 5.2851 +PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==) 5.2852 +PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=) 5.2853 +PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=) 5.2854 +PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >) 5.2855 +PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp) 5.2856 + 5.2857 + 5.2858 +/*--------------------------------* 5.2859 + * GiST index support functions * 5.2860 + *--------------------------------*/ 5.2861 + 5.2862 +/* GiST "consistent" support function */ 5.2863 +PG_FUNCTION_INFO_V1(pgl_gist_consistent); 5.2864 +Datum pgl_gist_consistent(PG_FUNCTION_ARGS) { 5.2865 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 5.2866 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 5.2867 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 5.2868 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 5.2869 + /* demand recheck because index and query methods are lossy */ 5.2870 + *recheck = true; 5.2871 + /* strategy number aliases for different operators using the same strategy */ 5.2872 + strategy %= 100; 5.2873 + /* strategy number 11: equality of two points */ 5.2874 + if (strategy == 11) { 5.2875 + /* query datum is another point */ 5.2876 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 5.2877 + /* convert other point to key */ 5.2878 + pgl_pointkey querykey; 5.2879 + pgl_point_to_key(query, querykey); 5.2880 + /* return true if both keys overlap */ 5.2881 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 5.2882 + } 5.2883 + /* strategy number 13: equality of two circles */ 5.2884 + if (strategy == 13) { 5.2885 + /* query datum is another circle */ 5.2886 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 5.2887 + /* convert other circle to key */ 5.2888 + pgl_areakey querykey; 5.2889 + pgl_circle_to_key(query, querykey); 5.2890 + /* return true if both keys overlap */ 5.2891 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 5.2892 + } 5.2893 + /* for all remaining strategies, keys on empty objects produce no match */ 5.2894 + /* (check necessary because query radius may be infinite) */ 5.2895 + if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false); 5.2896 + /* strategy number 21: overlapping with point */ 5.2897 + if (strategy == 21) { 5.2898 + /* query datum is a point */ 5.2899 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 5.2900 + /* return true if estimated distance (allowed to be smaller than real 5.2901 + distance) between index key and point is zero */ 5.2902 + PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0); 5.2903 + } 5.2904 + /* strategy number 22: (point) overlapping with box */ 5.2905 + if (strategy == 22) { 5.2906 + /* query datum is a box */ 5.2907 + pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1); 5.2908 + /* determine bounding box of indexed key */ 5.2909 + pgl_box keybox; 5.2910 + pgl_key_to_box(key, &keybox); 5.2911 + /* return true if query box overlaps with bounding box of indexed key */ 5.2912 + PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox)); 5.2913 + } 5.2914 + /* strategy number 23: overlapping with circle */ 5.2915 + if (strategy == 23) { 5.2916 + /* query datum is a circle */ 5.2917 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 5.2918 + /* return true if estimated distance (allowed to be smaller than real 5.2919 + distance) between index key and circle center is smaller than radius */ 5.2920 + PG_RETURN_BOOL( 5.2921 + pgl_estimate_key_distance(key, &(query->center)) <= query->radius 5.2922 + ); 5.2923 + } 5.2924 + /* strategy number 24: overlapping with cluster */ 5.2925 + if (strategy == 24) { 5.2926 + bool retval; /* return value */ 5.2927 + /* query datum is a cluster */ 5.2928 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.2929 + /* return true if estimated distance (allowed to be smaller than real 5.2930 + distance) between index key and circle center is smaller than radius */ 5.2931 + retval = ( 5.2932 + pgl_estimate_key_distance(key, &(query->bounding.center)) <= 5.2933 + query->bounding.radius 5.2934 + ); 5.2935 + PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */ 5.2936 + PG_RETURN_BOOL(retval); 5.2937 + } 5.2938 + /* throw error for any unknown strategy number */ 5.2939 + elog(ERROR, "unrecognized strategy number: %d", strategy); 5.2940 +} 5.2941 + 5.2942 +/* GiST "union" support function */ 5.2943 +PG_FUNCTION_INFO_V1(pgl_gist_union); 5.2944 +Datum pgl_gist_union(PG_FUNCTION_ARGS) { 5.2945 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 5.2946 + pgl_keyptr out; /* return value (to be palloc'ed) */ 5.2947 + int i; 5.2948 + /* determine key size */ 5.2949 + size_t keysize = PGL_KEY_IS_AREAKEY( 5.2950 + (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key) 5.2951 + ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey); 5.2952 + /* begin with first key as result */ 5.2953 + out = palloc(keysize); 5.2954 + memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize); 5.2955 + /* unite current result with second, third, etc. key */ 5.2956 + for (i=1; i<entryvec->n; i++) { 5.2957 + pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)); 5.2958 + } 5.2959 + /* return result */ 5.2960 + PG_RETURN_POINTER(out); 5.2961 +} 5.2962 + 5.2963 +/* GiST "compress" support function for indicis on points */ 5.2964 +PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint); 5.2965 +Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) { 5.2966 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 5.2967 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 5.2968 + /* only transform new leaves */ 5.2969 + if (entry->leafkey) { 5.2970 + /* get point to be transformed */ 5.2971 + pgl_point *point = (pgl_point *)DatumGetPointer(entry->key); 5.2972 + /* allocate memory for key */ 5.2973 + pgl_keyptr key = palloc(sizeof(pgl_pointkey)); 5.2974 + /* transform point to key */ 5.2975 + pgl_point_to_key(point, key); 5.2976 + /* create new GISTENTRY structure as return value */ 5.2977 + retval = palloc(sizeof(GISTENTRY)); 5.2978 + gistentryinit( 5.2979 + *retval, PointerGetDatum(key), 5.2980 + entry->rel, entry->page, entry->offset, FALSE 5.2981 + ); 5.2982 + } else { 5.2983 + /* inner nodes have already been transformed */ 5.2984 + retval = entry; 5.2985 + } 5.2986 + /* return pointer to old or new GISTENTRY structure */ 5.2987 + PG_RETURN_POINTER(retval); 5.2988 +} 5.2989 + 5.2990 +/* GiST "compress" support function for indicis on circles */ 5.2991 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle); 5.2992 +Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) { 5.2993 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 5.2994 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 5.2995 + /* only transform new leaves */ 5.2996 + if (entry->leafkey) { 5.2997 + /* get circle to be transformed */ 5.2998 + pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key); 5.2999 + /* allocate memory for key */ 5.3000 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 5.3001 + /* transform circle to key */ 5.3002 + pgl_circle_to_key(circle, key); 5.3003 + /* create new GISTENTRY structure as return value */ 5.3004 + retval = palloc(sizeof(GISTENTRY)); 5.3005 + gistentryinit( 5.3006 + *retval, PointerGetDatum(key), 5.3007 + entry->rel, entry->page, entry->offset, FALSE 5.3008 + ); 5.3009 + } else { 5.3010 + /* inner nodes have already been transformed */ 5.3011 + retval = entry; 5.3012 + } 5.3013 + /* return pointer to old or new GISTENTRY structure */ 5.3014 + PG_RETURN_POINTER(retval); 5.3015 +} 5.3016 + 5.3017 +/* GiST "compress" support function for indices on clusters */ 5.3018 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster); 5.3019 +Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) { 5.3020 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 5.3021 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 5.3022 + /* only transform new leaves */ 5.3023 + if (entry->leafkey) { 5.3024 + /* get cluster to be transformed (detoasting necessary!) */ 5.3025 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key); 5.3026 + /* allocate memory for key */ 5.3027 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 5.3028 + /* transform cluster to key */ 5.3029 + pgl_circle_to_key(&(cluster->bounding), key); 5.3030 + /* create new GISTENTRY structure as return value */ 5.3031 + retval = palloc(sizeof(GISTENTRY)); 5.3032 + gistentryinit( 5.3033 + *retval, PointerGetDatum(key), 5.3034 + entry->rel, entry->page, entry->offset, FALSE 5.3035 + ); 5.3036 + /* free detoasted datum */ 5.3037 + if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster); 5.3038 + } else { 5.3039 + /* inner nodes have already been transformed */ 5.3040 + retval = entry; 5.3041 + } 5.3042 + /* return pointer to old or new GISTENTRY structure */ 5.3043 + PG_RETURN_POINTER(retval); 5.3044 +} 5.3045 + 5.3046 +/* GiST "decompress" support function for indices */ 5.3047 +PG_FUNCTION_INFO_V1(pgl_gist_decompress); 5.3048 +Datum pgl_gist_decompress(PG_FUNCTION_ARGS) { 5.3049 + /* return passed pointer without transformation */ 5.3050 + PG_RETURN_POINTER(PG_GETARG_POINTER(0)); 5.3051 +} 5.3052 + 5.3053 +/* GiST "penalty" support function */ 5.3054 +PG_FUNCTION_INFO_V1(pgl_gist_penalty); 5.3055 +Datum pgl_gist_penalty(PG_FUNCTION_ARGS) { 5.3056 + GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0); 5.3057 + GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1); 5.3058 + float *penalty = (float *)PG_GETARG_POINTER(2); 5.3059 + /* get original key and key to insert */ 5.3060 + pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key); 5.3061 + pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key); 5.3062 + /* copy original key */ 5.3063 + union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key; 5.3064 + if (PGL_KEY_IS_AREAKEY(orig)) { 5.3065 + memcpy(union_key.areakey, orig, sizeof(union_key.areakey)); 5.3066 + } else { 5.3067 + memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey)); 5.3068 + } 5.3069 + /* calculate union of both keys */ 5.3070 + pgl_unite_keys((pgl_keyptr)&union_key, new); 5.3071 + /* penalty equal to reduction of key length (logarithm of added area) */ 5.3072 + /* (return value by setting referenced value and returning pointer) */ 5.3073 + *penalty = ( 5.3074 + PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key) 5.3075 + ); 5.3076 + PG_RETURN_POINTER(penalty); 5.3077 +} 5.3078 + 5.3079 +/* GiST "picksplit" support function */ 5.3080 +PG_FUNCTION_INFO_V1(pgl_gist_picksplit); 5.3081 +Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) { 5.3082 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 5.3083 + GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1); 5.3084 + OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */ 5.3085 + union { 5.3086 + pgl_pointkey pointkey; 5.3087 + pgl_areakey areakey; 5.3088 + } union_all; /* union of all keys (to be calculated from scratch) 5.3089 + (later cut in half) */ 5.3090 + int is_areakey = PGL_KEY_IS_AREAKEY( 5.3091 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key) 5.3092 + ); 5.3093 + int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey); 5.3094 + pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */ 5.3095 + pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */ 5.3096 + pgl_keyptr key; /* current key to be processed */ 5.3097 + /* allocate memory for array of left and right keys, set counts to zero */ 5.3098 + v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 5.3099 + v->spl_nleft = 0; 5.3100 + v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 5.3101 + v->spl_nright = 0; 5.3102 + /* calculate union of all keys from scratch */ 5.3103 + memcpy( 5.3104 + (pgl_keyptr)&union_all, 5.3105 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key), 5.3106 + keysize 5.3107 + ); 5.3108 + for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) { 5.3109 + pgl_unite_keys( 5.3110 + (pgl_keyptr)&union_all, 5.3111 + (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key) 5.3112 + ); 5.3113 + } 5.3114 + /* check if trivial split is necessary due to exhausted key length */ 5.3115 + /* (Note: keys for empty objects must have node depth set to maximum) */ 5.3116 + if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == ( 5.3117 + is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH 5.3118 + )) { 5.3119 + /* half of all keys go left */ 5.3120 + for ( 5.3121 + i=FirstOffsetNumber; 5.3122 + i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 5.3123 + i=OffsetNumberNext(i) 5.3124 + ) { 5.3125 + /* pointer to current key */ 5.3126 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 5.3127 + /* update unionL */ 5.3128 + /* check if key is first key that goes left */ 5.3129 + if (!v->spl_nleft) { 5.3130 + /* first key that goes left is just copied to unionL */ 5.3131 + memcpy(unionL, key, keysize); 5.3132 + } else { 5.3133 + /* unite current value and next key */ 5.3134 + pgl_unite_keys(unionL, key); 5.3135 + } 5.3136 + /* append offset number to list of keys that go left */ 5.3137 + v->spl_left[v->spl_nleft++] = i; 5.3138 + } 5.3139 + /* other half goes right */ 5.3140 + for ( 5.3141 + i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 5.3142 + i<entryvec->n; 5.3143 + i=OffsetNumberNext(i) 5.3144 + ) { 5.3145 + /* pointer to current key */ 5.3146 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 5.3147 + /* update unionR */ 5.3148 + /* check if key is first key that goes right */ 5.3149 + if (!v->spl_nright) { 5.3150 + /* first key that goes right is just copied to unionR */ 5.3151 + memcpy(unionR, key, keysize); 5.3152 + } else { 5.3153 + /* unite current value and next key */ 5.3154 + pgl_unite_keys(unionR, key); 5.3155 + } 5.3156 + /* append offset number to list of keys that go right */ 5.3157 + v->spl_right[v->spl_nright++] = i; 5.3158 + } 5.3159 + } 5.3160 + /* otherwise, a non-trivial split is possible */ 5.3161 + else { 5.3162 + /* cut covered area in half */ 5.3163 + /* (union_all then refers to area of keys that go left) */ 5.3164 + /* check if union of all keys covers empty and non-empty objects */ 5.3165 + if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) { 5.3166 + /* if yes, split into empty and non-empty objects */ 5.3167 + pgl_key_set_empty((pgl_keyptr)&union_all); 5.3168 + } else { 5.3169 + /* otherwise split by next bit */ 5.3170 + ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++; 5.3171 + /* NOTE: type bit conserved */ 5.3172 + } 5.3173 + /* determine for each key if it goes left or right */ 5.3174 + for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) { 5.3175 + /* pointer to current key */ 5.3176 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 5.3177 + /* keys within one half of the area go left */ 5.3178 + if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) { 5.3179 + /* update unionL */ 5.3180 + /* check if key is first key that goes left */ 5.3181 + if (!v->spl_nleft) { 5.3182 + /* first key that goes left is just copied to unionL */ 5.3183 + memcpy(unionL, key, keysize); 5.3184 + } else { 5.3185 + /* unite current value of unionL and processed key */ 5.3186 + pgl_unite_keys(unionL, key); 5.3187 + } 5.3188 + /* append offset number to list of keys that go left */ 5.3189 + v->spl_left[v->spl_nleft++] = i; 5.3190 + } 5.3191 + /* the other keys go right */ 5.3192 + else { 5.3193 + /* update unionR */ 5.3194 + /* check if key is first key that goes right */ 5.3195 + if (!v->spl_nright) { 5.3196 + /* first key that goes right is just copied to unionR */ 5.3197 + memcpy(unionR, key, keysize); 5.3198 + } else { 5.3199 + /* unite current value of unionR and processed key */ 5.3200 + pgl_unite_keys(unionR, key); 5.3201 + } 5.3202 + /* append offset number to list of keys that go right */ 5.3203 + v->spl_right[v->spl_nright++] = i; 5.3204 + } 5.3205 + } 5.3206 + } 5.3207 + /* store unions in return value */ 5.3208 + v->spl_ldatum = PointerGetDatum(unionL); 5.3209 + v->spl_rdatum = PointerGetDatum(unionR); 5.3210 + /* return all results */ 5.3211 + PG_RETURN_POINTER(v); 5.3212 +} 5.3213 + 5.3214 +/* GiST "same"/"equal" support function */ 5.3215 +PG_FUNCTION_INFO_V1(pgl_gist_same); 5.3216 +Datum pgl_gist_same(PG_FUNCTION_ARGS) { 5.3217 + pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0); 5.3218 + pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1); 5.3219 + bool *boolptr = (bool *)PG_GETARG_POINTER(2); 5.3220 + /* two keys are equal if they are binary equal */ 5.3221 + /* (return result by setting referenced boolean and returning pointer) */ 5.3222 + *boolptr = !memcmp( 5.3223 + key1, 5.3224 + key2, 5.3225 + PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey) 5.3226 + ); 5.3227 + PG_RETURN_POINTER(boolptr); 5.3228 +} 5.3229 + 5.3230 +/* GiST "distance" support function */ 5.3231 +PG_FUNCTION_INFO_V1(pgl_gist_distance); 5.3232 +Datum pgl_gist_distance(PG_FUNCTION_ARGS) { 5.3233 + GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0); 5.3234 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 5.3235 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 5.3236 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 5.3237 + double distance; /* return value */ 5.3238 + /* demand recheck because distance is just an estimation */ 5.3239 + /* (real distance may be bigger) */ 5.3240 + *recheck = true; 5.3241 + /* strategy number aliases for different operators using the same strategy */ 5.3242 + strategy %= 100; 5.3243 + /* strategy number 31: distance to point */ 5.3244 + if (strategy == 31) { 5.3245 + /* query datum is a point */ 5.3246 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 5.3247 + /* use pgl_estimate_pointkey_distance() function to compute result */ 5.3248 + distance = pgl_estimate_key_distance(key, query); 5.3249 + /* avoid infinity (reserved!) */ 5.3250 + if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 5.3251 + /* return result */ 5.3252 + PG_RETURN_FLOAT8(distance); 5.3253 + } 5.3254 + /* strategy number 33: distance to circle */ 5.3255 + if (strategy == 33) { 5.3256 + /* query datum is a circle */ 5.3257 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 5.3258 + /* estimate distance to circle center and substract circle radius */ 5.3259 + distance = ( 5.3260 + pgl_estimate_key_distance(key, &(query->center)) - query->radius 5.3261 + ); 5.3262 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 5.3263 + if (distance <= 0) distance = 0; 5.3264 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 5.3265 + /* return result */ 5.3266 + PG_RETURN_FLOAT8(distance); 5.3267 + } 5.3268 + /* strategy number 34: distance to cluster */ 5.3269 + if (strategy == 34) { 5.3270 + /* query datum is a cluster */ 5.3271 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 5.3272 + /* estimate distance to bounding center and substract bounding radius */ 5.3273 + distance = ( 5.3274 + pgl_estimate_key_distance(key, &(query->bounding.center)) - 5.3275 + query->bounding.radius 5.3276 + ); 5.3277 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 5.3278 + if (distance <= 0) distance = 0; 5.3279 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 5.3280 + /* free detoasted cluster (if copy) */ 5.3281 + PG_FREE_IF_COPY(query, 1); 5.3282 + /* return result */ 5.3283 + PG_RETURN_FLOAT8(distance); 5.3284 + } 5.3285 + /* throw error for any unknown strategy number */ 5.3286 + elog(ERROR, "unrecognized strategy number: %d", strategy); 5.3287 +} 5.3288 +