diff --git a/Makefile b/Makefile
index 5f428ea..1e9805f 100644
--- a/Makefile
+++ b/Makefile
@@ -10,8 +10,8 @@ SRC_DIR = $(shell basename $(shell pwd))
MODULE_big = pg_sphere
OBJS = src/sscan.o src/sparse.o src/sbuffer.o src/vector3d.o src/point.o \
- src/euler.o src/circle.o src/line.o src/ellipse.o src/polygon.o \
- src/path.o src/box.o src/output.o src/gq_cache.o src/gist.o \
+ src/euler.o src/circle.o src/circle_sel.o src/line.o src/ellipse.o src/polygon.o \
+ src/path.o src/box.o src/output.o src/gq_cache.o src/gist.o src/gist_support.o \
src/key.o src/gnomo.o src/epochprop.o src/brin.o
ifneq ($(USE_HEALPIX),0)
@@ -35,11 +35,11 @@ DATA_built = $(RELEASE_SQL) \
DOCS = README.pg_sphere COPYRIGHT.pg_sphere
REGRESS = init tables points euler circle line ellipse poly path box index \
contains_ops contains_ops_compat bounding_box_gist gnomo epochprop \
- contains overlaps spoint_brin sbox_brin
+ contains overlaps spoint_brin sbox_brin selectivity
TESTS = init_test tables points euler circle line ellipse poly path box \
index contains_ops contains_ops_compat bounding_box_gist gnomo \
- epochprop contains overlaps spoint_brin sbox_brin
+ epochprop contains overlaps spoint_brin sbox_brin selectivity
PG_CFLAGS += -DPGSPHERE_VERSION=$(PGSPHERE_VERSION)
PG_CPPFLAGS += -DPGSPHERE_VERSION=$(PGSPHERE_VERSION)
@@ -58,7 +58,7 @@ CRUSH_TESTS = init_extended circle_extended
PGS_SQL = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \
pgs_line.sql pgs_ellipse.sql pgs_polygon.sql pgs_path.sql \
pgs_box.sql pgs_contains_ops.sql pgs_contains_ops_compat.sql \
- pgs_gist.sql gnomo.sql pgs_brin.sql
+ pgs_gist.sql gnomo.sql pgs_brin.sql pgs_circle_sel.sql
ifneq ($(USE_HEALPIX),0)
REGRESS += healpix moc moc1 moc100 mocautocast
@@ -102,10 +102,17 @@ healpix_bare/healpix_bare.o : healpix_bare/healpix_bare.c
$(COMPILE.c) -Wno-declaration-after-statement -o $@ $^
pg_version := $(word 2,$(shell $(PG_CONFIG) --version))
+has_support_functions = $(if $(filter-out 9.% 10.% 11.%,$(pg_version)),y,n)
crushtest: REGRESS += $(CRUSH_TESTS)
crushtest: installcheck
+ifeq ($(has_support_functions),y)
+PGS_SQL += pgs_gist_support.sql
+REGRESS += gist_support
+TESTS += gist_support
+endif
+
test: pg_sphere.test.sql
$(pg_regress_installcheck) --temp-instance=tmp_check $(REGRESS_OPTS) $(TESTS)
@@ -180,8 +187,11 @@ pg_sphere--1.2.3--1.3.0.sql: pgs_brin.sql.in
pg_sphere--1.3.0--1.3.1.sql:
cat upgrade_scripts/$@.in > $@
-pg_sphere--1.3.1--1.3.2.sql:
- cat upgrade_scripts/$@.in > $@
+ifeq ($(has_support_functions),y)
+pg_sphere--1.3.1--1.3.2.sql: pgs_gist_support.sql.in
+endif
+pg_sphere--1.3.1--1.3.2.sql: pgs_circle_sel.sql.in
+ cat upgrade_scripts/$@.in $^ > $@
# end of local stuff
diff --git a/doc/functions.sgm b/doc/functions.sgm
index b65f8c0..232b755 100644
--- a/doc/functions.sgm
+++ b/doc/functions.sgm
@@ -149,6 +149,41 @@
+
+
+ Point-within-distance function
+
+
+ The function
+
+
+
+ spoint_dwithin
+ spoint p1
+ spoint p2
+ float8 radius
+
+
+
+ returns a boolean value that signifies whether the points
+ p1 and p2
+ lie within distance radius (in radians) of each other, i.e.
+ it computes the boolean expression p1 <-> p2 <= radius.
+ On PostgreSQL 12 and later, the function has GiST
+ support and the PostgreSQL optimizer will transform it to either
+ p1 <@ scircle(p2, radius) or
+ p2 <@ scircle(p1, radius) where appropriate.
+
+
+
+ Efficiently join two tables of points with some fuzziness permitted
+
+
+ SELECT * FROM stars1 JOIN stars2 WHERE spoint_dwithin(stars1.s, stars2.s, 1e-5);]]>
+
+
+
+
diff --git a/expected/gist_support.out b/expected/gist_support.out
new file mode 100644
index 0000000..dead655
--- /dev/null
+++ b/expected/gist_support.out
@@ -0,0 +1,155 @@
+-- spoint_dwithin function selectivity
+set jit = off; -- suppress extra planning output
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), 1)');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '1'::double precision)
+ Rows Removed by Filter: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), .1)');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '0.1'::double precision)
+ Rows Removed by Filter: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), .01)');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.01>'::scircle)
+(2 rows)
+
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, 1)');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Filter: spoint_dwithin('(1 , 1)'::spoint, star, '1'::double precision)
+ Rows Removed by Filter: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, .1)');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Filter: spoint_dwithin('(1 , 1)'::spoint, star, '0.1'::double precision)
+ Rows Removed by Filter: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, .01)');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.01>'::scircle)
+(2 rows)
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, 1)', do_analyze := 'false');
+ explain
+---------------------------------------------------------------------------------------
+ Nested Loop (rows=22984885 width=32)
+ -> Seq Scan on spoint10k a (rows=10000 width=16)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=2298 width=16)
+ Index Cond: (star OPERATOR(public.<@) scircle(a.star, '1'::double precision))
+(4 rows)
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .1)');
+ explain
+-----------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=249792 width=32) (actual rows=505342 loops=1)
+ -> Seq Scan on spoint10k a (rows=10000 width=16) (actual rows=10000 loops=1)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=25 width=16) (actual rows=51 loops=10000)
+ Index Cond: (star OPERATOR(public.<@) scircle(a.star, '0.1'::double precision))
+ Rows Removed by Index Recheck: 31
+(5 rows)
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)');
+ explain
+---------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=2500 width=32) (actual rows=17614 loops=1)
+ -> Seq Scan on spoint10k a (rows=10000 width=16) (actual rows=10000 loops=1)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=1 width=16) (actual rows=2 loops=10000)
+ Index Cond: (star OPERATOR(public.<@) scircle(a.star, '0.01'::double precision))
+ Rows Removed by Index Recheck: 1
+(5 rows)
+
+-- spoint_dwithin is symmetric in the first two arguments
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1)');
+ explain
+------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=6 width=32) (actual rows=33 loops=1)
+ -> Bitmap Heap Scan on spoint10k a (rows=25 width=16) (actual rows=29 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '0.1'::double precision)
+ Rows Removed by Filter: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=1 width=16) (actual rows=1 loops=29)
+ Index Cond: (star OPERATOR(public.<@) scircle(a.star, '0.01'::double precision))
+ Rows Removed by Index Recheck: 0
+(10 rows)
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(b.star, a.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1)');
+ explain
+------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=6 width=32) (actual rows=33 loops=1)
+ -> Bitmap Heap Scan on spoint10k a (rows=25 width=16) (actual rows=29 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '0.1'::double precision)
+ Rows Removed by Filter: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=1 width=16) (actual rows=1 loops=29)
+ Index Cond: (star OPERATOR(public.<@) scircle(a.star, '0.01'::double precision))
+ Rows Removed by Index Recheck: 0
+(10 rows)
+
+-- both sides indexable, check if the planner figures out the better choice
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1) and spoint_dwithin(b.star, spoint(1,1), .05)');
+ explain
+-------------------------------------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=1 width=32) (actual rows=16 loops=1)
+ -> Bitmap Heap Scan on spoint10k b (rows=6 width=16) (actual rows=12 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '0.05'::double precision)
+ Rows Removed by Filter: 4
+ Heap Blocks: exact=14
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=6 width=0) (actual rows=16 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.05>'::scircle)
+ -> Index Scan using spoint10k_star_idx on spoint10k a (rows=1 width=16) (actual rows=1 loops=12)
+ Index Cond: ((star OPERATOR(public.<@) scircle(b.star, '0.01'::double precision)) AND (star <@ '<(1 , 1) , 0.1>'::scircle))
+ Rows Removed by Index Recheck: 0
+(10 rows)
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .05) and spoint_dwithin(b.star, spoint(1,1), .1)');
+ explain
+-------------------------------------------------------------------------------------------------------------------------------------
+ Nested Loop (rows=1 width=32) (actual rows=16 loops=1)
+ -> Bitmap Heap Scan on spoint10k a (rows=6 width=16) (actual rows=12 loops=1)
+ Filter: spoint_dwithin(star, '(1 , 1)'::spoint, '0.05'::double precision)
+ Rows Removed by Filter: 4
+ Heap Blocks: exact=14
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=6 width=0) (actual rows=16 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.05>'::scircle)
+ -> Index Scan using spoint10k_star_idx on spoint10k b (rows=1 width=16) (actual rows=1 loops=12)
+ Index Cond: ((star OPERATOR(public.<@) scircle(a.star, '0.01'::double precision)) AND (star <@ '<(1 , 1) , 0.1>'::scircle))
+ Rows Removed by Index Recheck: 0
+(10 rows)
+
diff --git a/expected/index.out b/expected/index.out
index 8c67155..eefa5c6 100644
--- a/expected/index.out
+++ b/expected/index.out
@@ -55,6 +55,7 @@ SELECT count(*) FROM spheretmp4 WHERE l && scircle '<(1,1),0.3>';
-- create idx
CREATE TABLE spheretmp1b AS TABLE spheretmp1;
+ANALYZE spheretmp1;
CREATE INDEX aaaidx ON spheretmp1 USING gist ( p );
CREATE INDEX spoint3_idx ON spheretmp1b USING gist (p spoint3);
CREATE INDEX bbbidx ON spheretmp2 USING gist ( c );
@@ -165,7 +166,6 @@ EXPLAIN (COSTS OFF) SELECT count(*) FROM spheretmp1b WHERE p = spoint '(3.09 , 1
4
(1 row)
-SET enable_bitmapscan = ON;
SET enable_indexonlyscan = OFF;
EXPLAIN (COSTS OFF) SELECT count(*) FROM spheretmp1b WHERE p <@ scircle '<(1,1),0.3>';
QUERY PLAN
diff --git a/expected/points.out b/expected/points.out
index 6fb3c17..f3a0713 100644
--- a/expected/points.out
+++ b/expected/points.out
@@ -666,3 +666,19 @@ SELECT '( 0h 2m 30s , -90d 0m 0s)'::spoint<->'( 12h 2m 30s , -90d 0m 0s)'::spoin
0
(1 row)
+-- spoint_dwithin function ----------
+SELECT a, b, radius, a <-> b AS "<->", spoint_dwithin(a, b, radius)
+FROM (VALUES
+ ('(0, 0)'::spoint, '(0, 0)'::spoint, 0),
+ ('(0, 0)', '(0, 1)', 1),
+ ('(0, 0)', '(0.1, 0.1)', 0.14),
+ ('(0, 0)', '(0.1, 0.1)', 0.15)
+ ) sub (a, b, radius);
+ a | b | radius | <-> | spoint_dwithin
+---------+-------------+--------+-----------------+----------------
+ (0 , 0) | (0 , 0) | 0 | 0 | t
+ (0 , 0) | (0 , 1) | 1 | 1 | t
+ (0 , 0) | (0.1 , 0.1) | 0.14 | 0.1413032986961 | f
+ (0 , 0) | (0.1 , 0.1) | 0.15 | 0.1413032986961 | t
+(4 rows)
+
diff --git a/expected/selectivity.out b/expected/selectivity.out
new file mode 100644
index 0000000..e694f0f
--- /dev/null
+++ b/expected/selectivity.out
@@ -0,0 +1,124 @@
+-- test selectivity estimator functions
+create table spoint10k (star spoint);
+insert into spoint10k select spoint(i, i*i) from generate_series(1, 10000) g(i);
+create index on spoint10k using gist (star);
+analyze spoint10k;
+-- "explain analyze" wrapper that removes 'cost=...' since it varies across architectures
+-- (we can't use "costs off" since that also removes the estimated row count)
+create or replace function explain(query text, do_analyze text default 'true') returns setof text language plpgsql as $$
+declare
+ line text;
+begin
+ for line in execute format('explain (analyze %s, timing off, summary off) %s', do_analyze, query) loop
+ return next regexp_replace(line, 'cost=\S+ ', '');
+ end loop;
+ return;
+end;
+$$;
+-- <@ operator selectivity
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), 1)');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Recheck Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+ Rows Removed by Index Recheck: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .1)');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Recheck Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+ Rows Removed by Index Recheck: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .01)');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.01>'::scircle)
+(2 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) @> star');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Recheck Cond: ('<(1 , 1) , 1>'::scircle @> star)
+ Rows Removed by Index Recheck: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) @> star');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Recheck Cond: ('<(1 , 1) , 0.1>'::scircle @> star)
+ Rows Removed by Index Recheck: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) @> star');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.01>'::scircle)
+(2 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), 1)');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=7702 width=16) (actual rows=6991 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 1>'::scircle)
+ Rows Removed by Filter: 3009
+(3 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .1)');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=9975 width=16) (actual rows=9971 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 0.1>'::scircle)
+ Rows Removed by Filter: 29
+(3 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .01)');
+ explain
+-------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=10000 width=16) (actual rows=9999 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 0.01>'::scircle)
+ Rows Removed by Filter: 1
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) !@> star');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=7702 width=16) (actual rows=6991 loops=1)
+ Filter: ('<(1 , 1) , 1>'::scircle !@> star)
+ Rows Removed by Filter: 3009
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) !@> star');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=9975 width=16) (actual rows=9971 loops=1)
+ Filter: ('<(1 , 1) , 0.1>'::scircle !@> star)
+ Rows Removed by Filter: 29
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) !@> star');
+ explain
+-------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=10000 width=16) (actual rows=9999 loops=1)
+ Filter: ('<(1 , 1) , 0.01>'::scircle !@> star)
+ Rows Removed by Filter: 1
+(3 rows)
+
diff --git a/expected/selectivity_1.out b/expected/selectivity_1.out
new file mode 100644
index 0000000..08cd8f7
--- /dev/null
+++ b/expected/selectivity_1.out
@@ -0,0 +1,124 @@
+-- test selectivity estimator functions
+create table spoint10k (star spoint);
+insert into spoint10k select spoint(i, i*i) from generate_series(1, 10000) g(i);
+create index on spoint10k using gist (star);
+analyze spoint10k;
+-- "explain analyze" wrapper that removes 'cost=...' since it varies across architectures
+-- (we can't use "costs off" since that also removes the estimated row count)
+create or replace function explain(query text, do_analyze text default 'true') returns setof text language plpgsql as $$
+declare
+ line text;
+begin
+ for line in execute format('explain (analyze %s, timing off, summary off) %s', do_analyze, query) loop
+ return next regexp_replace(line, 'cost=\S+ ', '');
+ end loop;
+ return;
+end;
+$$;
+-- <@ operator selectivity
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), 1)');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Recheck Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+ Rows Removed by Index Recheck: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .1)');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Recheck Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+ Rows Removed by Index Recheck: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.1>'::scircle)
+(6 rows)
+
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .01)');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: (star <@ '<(1 , 1) , 0.01>'::scircle)
+(2 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) @> star');
+ explain
+-----------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=2298 width=16) (actual rows=3009 loops=1)
+ Recheck Cond: ('<(1 , 1) , 1>'::scircle @> star)
+ Rows Removed by Index Recheck: 1560
+ Heap Blocks: exact=55
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=2298 width=0) (actual rows=4569 loops=1)
+ Index Cond: ('<(1 , 1) , 1>'::scircle @> star)
+(6 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) @> star');
+ explain
+-------------------------------------------------------------------------------------------
+ Bitmap Heap Scan on spoint10k (rows=25 width=16) (actual rows=29 loops=1)
+ Recheck Cond: ('<(1 , 1) , 0.1>'::scircle @> star)
+ Rows Removed by Index Recheck: 19
+ Heap Blocks: exact=32
+ -> Bitmap Index Scan on spoint10k_star_idx (rows=25 width=0) (actual rows=48 loops=1)
+ Index Cond: ('<(1 , 1) , 0.1>'::scircle @> star)
+(6 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) @> star');
+ explain
+---------------------------------------------------------------------------------------------
+ Index Scan using spoint10k_star_idx on spoint10k (rows=1 width=16) (actual rows=1 loops=1)
+ Index Cond: ('<(1 , 1) , 0.01>'::scircle @> star)
+(2 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), 1)');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=7702 width=16) (actual rows=6991 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 1>'::scircle)
+ Rows Removed by Filter: 3009
+(3 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .1)');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=9975 width=16) (actual rows=9971 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 0.1>'::scircle)
+ Rows Removed by Filter: 29
+(3 rows)
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .01)');
+ explain
+-------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=10000 width=16) (actual rows=9999 loops=1)
+ Filter: (star !<@ '<(1 , 1) , 0.01>'::scircle)
+ Rows Removed by Filter: 1
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) !@> star');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=7702 width=16) (actual rows=6991 loops=1)
+ Filter: ('<(1 , 1) , 1>'::scircle !@> star)
+ Rows Removed by Filter: 3009
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) !@> star');
+ explain
+------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=9975 width=16) (actual rows=9971 loops=1)
+ Filter: ('<(1 , 1) , 0.1>'::scircle !@> star)
+ Rows Removed by Filter: 29
+(3 rows)
+
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) !@> star');
+ explain
+-------------------------------------------------------------------------
+ Seq Scan on spoint10k (rows=10000 width=16) (actual rows=9999 loops=1)
+ Filter: ('<(1 , 1) , 0.01>'::scircle !@> star)
+ Rows Removed by Filter: 1
+(3 rows)
+
diff --git a/pgs_circle_sel.sql.in b/pgs_circle_sel.sql.in
new file mode 100644
index 0000000..d1e9c7e
--- /dev/null
+++ b/pgs_circle_sel.sql.in
@@ -0,0 +1,109 @@
+CREATE FUNCTION spoint_contained_by_circle_sel(internal, oid, internal, integer)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_sel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_sel(internal, oid, internal, integer) IS
+ 'selectivity estimator function for spoint_contained_by_circle';
+
+CREATE FUNCTION spoint_contained_by_circle_joinsel(internal, oid, internal, smallint, internal)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_joinsel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_joinsel(internal, oid, internal, smallint, internal) IS
+ 'join selectivity estimator function for spoint_contained_by_circle';
+
+
+CREATE FUNCTION spoint_contained_by_circle_neg_sel(internal, oid, internal, integer)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_neg_sel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_neg_sel(internal, oid, internal, integer) IS
+ 'selectivity estimator function for spoint_contained_by_circle_neg';
+
+CREATE FUNCTION spoint_contained_by_circle_neg_joinsel(internal, oid, internal, smallint, internal)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_neg_joinsel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_neg_joinsel(internal, oid, internal, smallint, internal) IS
+ 'join selectivity estimator function for spoint_contained_by_circle_neg';
+
+
+CREATE FUNCTION spoint_contained_by_circle_com_sel(internal, oid, internal, integer)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_com_sel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_com_sel(internal, oid, internal, integer) IS
+ 'selectivity estimator function for spoint_contained_by_circle_com';
+
+CREATE FUNCTION spoint_contained_by_circle_com_joinsel(internal, oid, internal, smallint, internal)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_com_joinsel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_com_joinsel(internal, oid, internal, smallint, internal) IS
+ 'join selectivity estimator function for spoint_contained_by_circle_com';
+
+
+CREATE FUNCTION spoint_contained_by_circle_com_neg_sel(internal, oid, internal, integer)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_com_neg_sel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_com_neg_sel(internal, oid, internal, integer) IS
+ 'selectivity estimator function for spoint_contained_by_circle_com_neg';
+
+CREATE FUNCTION spoint_contained_by_circle_com_neg_joinsel(internal, oid, internal, smallint, internal)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME' , 'spherepoint_in_circle_com_neg_joinsel'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_contained_by_circle_com_neg_joinsel(internal, oid, internal, smallint, internal) IS
+ 'join selectivity estimator function for spoint_contained_by_circle_com_neg';
+
+
+ALTER OPERATOR @> (scircle, spoint)
+SET (
+ RESTRICT = spoint_contained_by_circle_com_sel,
+ JOIN = spoint_contained_by_circle_com_joinsel
+);
+
+ALTER OPERATOR <@ (spoint, scircle)
+SET (
+ RESTRICT = spoint_contained_by_circle_sel,
+ JOIN = spoint_contained_by_circle_joinsel
+);
+
+ALTER OPERATOR !@> (scircle, spoint)
+SET (
+ RESTRICT = spoint_contained_by_circle_com_neg_sel,
+ JOIN = spoint_contained_by_circle_com_neg_joinsel
+);
+
+ALTER OPERATOR !<@ (spoint, scircle)
+SET (
+ RESTRICT = spoint_contained_by_circle_neg_sel,
+ JOIN = spoint_contained_by_circle_neg_joinsel
+);
+
+
+CREATE FUNCTION spoint_dwithin(spoint, spoint, float8)
+ RETURNS boolean
+ AS 'MODULE_PATHNAME', 'spherepoint_dwithin'
+ LANGUAGE C
+ IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION spoint_dwithin(spoint, spoint, float8) IS
+ 'predicate matching spherical points less than a given distance apart';
diff --git a/pgs_gist_support.sql.in b/pgs_gist_support.sql.in
new file mode 100644
index 0000000..97b4059
--- /dev/null
+++ b/pgs_gist_support.sql.in
@@ -0,0 +1,12 @@
+-- PG12+ has support functions
+
+CREATE FUNCTION spoint_dwithin_supportfn (internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME', 'spherepoint_dwithin_supportfn'
+ LANGUAGE 'c';
+
+COMMENT ON FUNCTION spoint_dwithin_supportfn(internal) IS
+ 'support function for spoint_dwithin';
+
+ALTER FUNCTION spoint_dwithin(spoint, spoint, float8)
+ SUPPORT spoint_dwithin_supportfn;
diff --git a/sql/gist_support.sql b/sql/gist_support.sql
new file mode 100644
index 0000000..361b3d6
--- /dev/null
+++ b/sql/gist_support.sql
@@ -0,0 +1,26 @@
+-- spoint_dwithin function selectivity
+set jit = off; -- suppress extra planning output
+
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), 1)');
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), .1)');
+select explain('select * from spoint10k where spoint_dwithin(star, spoint(1,1), .01)');
+
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, 1)');
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, .1)');
+select explain('select * from spoint10k where spoint_dwithin(spoint(1,1), star, .01)');
+
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, 1)', do_analyze := 'false');
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .1)');
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)');
+
+-- spoint_dwithin is symmetric in the first two arguments
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1)');
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(b.star, a.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1)');
+
+-- both sides indexable, check if the planner figures out the better choice
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .1) and spoint_dwithin(b.star, spoint(1,1), .05)');
+select explain('select * from spoint10k a join spoint10k b on spoint_dwithin(a.star, b.star, .01)
+ where spoint_dwithin(a.star, spoint(1,1), .05) and spoint_dwithin(b.star, spoint(1,1), .1)');
diff --git a/sql/index.sql b/sql/index.sql
index fb3730a..8713106 100644
--- a/sql/index.sql
+++ b/sql/index.sql
@@ -24,6 +24,7 @@ SELECT count(*) FROM spheretmp4 WHERE l && scircle '<(1,1),0.3>';
-- create idx
CREATE TABLE spheretmp1b AS TABLE spheretmp1;
+ANALYZE spheretmp1;
CREATE INDEX aaaidx ON spheretmp1 USING gist ( p );
@@ -69,7 +70,6 @@ EXPLAIN (COSTS OFF) SELECT count(*) FROM spheretmp1b WHERE p <@ scircle '<(1,1),
EXPLAIN (COSTS OFF) SELECT count(*) FROM spheretmp1b WHERE p = spoint '(3.09 , 1.25)';
SELECT count(*) FROM spheretmp1b WHERE p = spoint '(3.09 , 1.25)';
-SET enable_bitmapscan = ON;
SET enable_indexonlyscan = OFF;
EXPLAIN (COSTS OFF) SELECT count(*) FROM spheretmp1b WHERE p <@ scircle '<(1,1),0.3>';
diff --git a/sql/points.sql b/sql/points.sql
index 5e63260..b34ae6e 100644
--- a/sql/points.sql
+++ b/sql/points.sql
@@ -240,3 +240,12 @@ SELECT '( 0h 2m 30s , 90d 0m 0s)'::spoint<->'( 12h 2m 30s , 90d 0m 0s)'::spoint;
SELECT '( 0h 2m 30s , -90d 0m 0s)'::spoint<->'( 12h 2m 30s , -90d 0m 0s)'::spoint;
+-- spoint_dwithin function ----------
+
+SELECT a, b, radius, a <-> b AS "<->", spoint_dwithin(a, b, radius)
+FROM (VALUES
+ ('(0, 0)'::spoint, '(0, 0)'::spoint, 0),
+ ('(0, 0)', '(0, 1)', 1),
+ ('(0, 0)', '(0.1, 0.1)', 0.14),
+ ('(0, 0)', '(0.1, 0.1)', 0.15)
+ ) sub (a, b, radius);
diff --git a/sql/selectivity.sql b/sql/selectivity.sql
new file mode 100644
index 0000000..d0f1d6e
--- /dev/null
+++ b/sql/selectivity.sql
@@ -0,0 +1,36 @@
+-- test selectivity estimator functions
+
+create table spoint10k (star spoint);
+insert into spoint10k select spoint(i, i*i) from generate_series(1, 10000) g(i);
+create index on spoint10k using gist (star);
+analyze spoint10k;
+
+-- "explain analyze" wrapper that removes 'cost=...' since it varies across architectures
+-- (we can't use "costs off" since that also removes the estimated row count)
+create or replace function explain(query text, do_analyze text default 'true') returns setof text language plpgsql as $$
+declare
+ line text;
+begin
+ for line in execute format('explain (analyze %s, timing off, summary off) %s', do_analyze, query) loop
+ return next regexp_replace(line, 'cost=\S+ ', '');
+ end loop;
+ return;
+end;
+$$;
+
+-- <@ operator selectivity
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), 1)');
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .1)');
+select explain('select * from spoint10k where star <@ scircle(spoint(1,1), .01)');
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) @> star');
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) @> star');
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) @> star');
+
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), 1)');
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .1)');
+select explain('select * from spoint10k where star !<@ scircle(spoint(1,1), .01)');
+
+select explain('select * from spoint10k where scircle(spoint(1,1), 1) !@> star');
+select explain('select * from spoint10k where scircle(spoint(1,1), .1) !@> star');
+select explain('select * from spoint10k where scircle(spoint(1,1), .01) !@> star');
diff --git a/src/circle.c b/src/circle.c
index e28fe0c..577ac3c 100644
--- a/src/circle.c
+++ b/src/circle.c
@@ -391,7 +391,7 @@ spherecircle_area(PG_FUNCTION_ARGS)
{
SCIRCLE *c = (SCIRCLE *) PG_GETARG_POINTER(0);
- PG_RETURN_FLOAT8(PID * (1 - cos(c->radius)));
+ PG_RETURN_FLOAT8(spherecircle_area_float(c->radius));
}
Datum
diff --git a/src/circle_sel.c b/src/circle_sel.c
new file mode 100644
index 0000000..d7bd843
--- /dev/null
+++ b/src/circle_sel.c
@@ -0,0 +1,218 @@
+#include "circle.h"
+#include
+#include
+
+/* Circle selectivity functions */
+
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_sel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_joinsel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_neg_sel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_neg_joinsel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_com_sel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_com_joinsel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_com_neg_sel);
+PG_FUNCTION_INFO_V1(spherepoint_in_circle_com_neg_joinsel);
+
+/*
+ * Common code for spherepoint_in_circle_sel()
+ */
+static double
+spherepoint_in_circle_sel_funcexpr(Node *other)
+{
+ FuncExpr *ofunc = (FuncExpr *)other;
+ char *func_name = get_func_name(ofunc->funcid);
+ Const *arg2;
+ double radius;
+ double selec;
+
+ /*
+ * Safety checks: are we called on a function called scircle that takes a
+ * const double as 2nd argument?
+ */
+ if (strcmp(func_name, "scircle") != 0)
+ {
+ ereport(DEBUG1, (errmsg("<@ selectivity called on function that is not scircle")));
+ return DEFAULT_SCIRCLE_SEL;
+ }
+ if (list_length(ofunc->args) != 2)
+ {
+ ereport(DEBUG1, (errmsg("<@ selectivity called on function that does not have 2 arguments")));
+ return DEFAULT_SCIRCLE_SEL;
+ }
+ if (!IsA(lsecond(ofunc->args), Const))
+ {
+ ereport(DEBUG1, (errmsg("<@ selectivity called on scircle() with non-const 2nd arg")));
+ return DEFAULT_SCIRCLE_SEL;
+ }
+ arg2 = (Const *) lsecond(ofunc->args);
+ if (arg2->consttype != FLOAT8OID)
+ {
+ ereport(DEBUG1, (errmsg("<@ selectivity called on scircle() with non-float8 2nd arg")));
+ return DEFAULT_SCIRCLE_SEL;
+ }
+
+ radius = DatumGetFloat8(arg2->constvalue);
+ selec = spherecircle_area_float(radius) / SPHERE_SURFACE;
+ CLAMP_PROBABILITY(selec);
+
+ return selec;
+}
+
+static double
+spherepoint_in_circle_sel_internal(PG_FUNCTION_ARGS, bool commute, bool negate)
+{
+ PlannerInfo *root = (PlannerInfo *) PG_GETARG_POINTER(0);
+ Oid operator = PG_GETARG_OID(1);
+ List *args = (List *) PG_GETARG_POINTER(2);
+ int varRelid = PG_GETARG_INT32(3);
+ //Oid collation = PG_GET_COLLATION();
+ VariableStatData vardata;
+ Node *other;
+ bool varonleft;
+ double selec;
+
+ /*
+ * When asked about <>, we do the estimation using the corresponding =
+ * operator, then convert to <> via "1.0 - eq_selectivity - nullfrac".
+ */
+ if (negate)
+ {
+ operator = get_negator(operator);
+ if (!OidIsValid(operator))
+ {
+ /* Use default selectivity (should we raise an error instead?) */
+ return 1.0 - DEFAULT_SCIRCLE_SEL;
+ }
+ }
+
+ /*
+ * If expression is not variable = something or something = variable, then
+ * punt and return a default estimate.
+ */
+ if (!get_restriction_variable(root, args, varRelid,
+ &vardata, &other, &varonleft))
+ return negate ? (1.0 - DEFAULT_SCIRCLE_SEL) : DEFAULT_SCIRCLE_SEL;
+
+ /*
+ * We can do a lot better if the something is a constant. (Note: the
+ * Const might result from estimation rather than being a simple constant
+ * in the query.)
+ * Look only at var op circle_const, not var op point_const.
+ */
+ if (IsA(other, Const) && varonleft != commute)
+ {
+ Const *constnode = (Const *) other;
+ SCIRCLE *c;
+
+ Assert(!constnode->constisnull); /* operators are strict, shouldn't have NULLs here */
+ c = (SCIRCLE *) constnode->constvalue;
+ selec = spherecircle_area_float(c->radius) / SPHERE_SURFACE;
+ CLAMP_PROBABILITY(selec);
+ }
+ /*
+ * Check for <@ scircle(..., radius)
+ */
+ else if (IsA(other, FuncExpr) && varonleft != commute)
+ {
+ selec = spherepoint_in_circle_sel_funcexpr(other);
+ // p *((double * )&((*(Const *)((((FuncExpr*) other)->args->elements)[1].ptr_value)).constvalue))
+ }
+ else
+ {
+ ereport(DEBUG1, (errmsg("<@ selectivity not const, other node tag %i", other->type)));
+ selec = DEFAULT_SCIRCLE_SEL;
+ }
+
+ ReleaseVariableStats(vardata);
+
+ return negate ? (1.0 - selec) : selec;
+}
+
+/*
+ * Common code for spherepoint_in_circle_joinsel()
+ */
+static double
+spherepoint_in_circle_joinsel_internal(PG_FUNCTION_ARGS, bool commute)
+{
+#ifdef NOT_USED
+ PlannerInfo *root = (PlannerInfo *) PG_GETARG_POINTER(0);
+ Oid operator = PG_GETARG_OID(1);
+ List *args = (List *) PG_GETARG_POINTER(2);
+
+ JoinType jointype = (JoinType) PG_GETARG_INT16(3);
+ SpecialJoinInfo *sjinfo = (SpecialJoinInfo *) PG_GETARG_POINTER(4);
+ Oid collation = PG_GET_COLLATION();
+ double selec;
+ double selec_inner;
+ VariableStatData vardata1;
+ VariableStatData vardata2;
+ double nd1;
+ double nd2;
+ bool isdefault1;
+ bool isdefault2;
+ Oid opfuncoid;
+ AttStatsSlot sslot1;
+ AttStatsSlot sslot2;
+ bool join_is_reversed;
+ RelOptInfo *inner_rel;
+
+ get_join_variables(root, args, sjinfo,
+ &vardata1, &vardata2, &join_is_reversed);
+
+ ReleaseVariableStats(vardata1);
+ ReleaseVariableStats(vardata2);
+#endif
+
+ /* return a constant default for now; practically joins should use the
+ * spoint_dwithin function instead which has join support with the
+ * additional advantage that it's symmetric */
+ return DEFAULT_SCIRCLE_SEL;
+}
+
+Datum
+spherepoint_in_circle_sel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_sel_internal(fcinfo, false, false));
+}
+
+Datum
+spherepoint_in_circle_joinsel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_joinsel_internal(fcinfo, false));
+}
+
+Datum
+spherepoint_in_circle_neg_sel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_sel_internal(fcinfo, false, true));
+}
+
+Datum
+spherepoint_in_circle_neg_joinsel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(DEFAULT_INEQ_SEL);
+}
+
+Datum
+spherepoint_in_circle_com_sel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_sel_internal(fcinfo, true, false));
+}
+
+Datum
+spherepoint_in_circle_com_joinsel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_joinsel_internal(fcinfo, true));
+}
+
+Datum
+spherepoint_in_circle_com_neg_sel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(spherepoint_in_circle_sel_internal(fcinfo, true, true));
+}
+
+Datum
+spherepoint_in_circle_com_neg_joinsel(PG_FUNCTION_ARGS)
+{
+ PG_RETURN_FLOAT8(DEFAULT_INEQ_SEL);
+}
diff --git a/src/gist_support.c b/src/gist_support.c
new file mode 100644
index 0000000..95aa65c
--- /dev/null
+++ b/src/gist_support.c
@@ -0,0 +1,223 @@
+/**********************************************************************
+ *
+ * pgsphere gist_support.c
+ * based on gserialized_supportfn.c from PostGIS
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS. If not, see .
+ *
+ **********************************************************************/
+
+
+/* PostgreSQL */
+#include "postgres.h"
+#if PG_VERSION_NUM >= 120000
+#include "funcapi.h"
+#include "access/htup_details.h"
+#include "access/stratnum.h"
+#include "catalog/namespace.h"
+#include "catalog/pg_opfamily.h"
+#include "catalog/pg_type_d.h"
+#include "catalog/pg_am_d.h"
+#include "nodes/supportnodes.h"
+#include "nodes/nodeFuncs.h"
+#include "nodes/makefuncs.h"
+#include "optimizer/optimizer.h"
+#include "parser/parse_func.h"
+#include "parser/parse_type.h"
+#include "utils/array.h"
+#include "utils/builtins.h"
+#include "utils/lsyscache.h"
+#include "utils/numeric.h"
+#include "utils/selfuncs.h"
+#include "utils/syscache.h"
+
+#include "point.h"
+#include "circle.h"
+
+static Oid
+scircleTypeOid(Oid callingfunc)
+{
+ /* type must be in same namespace as the caller */
+ char *nspname = get_namespace_name(get_func_namespace(callingfunc));
+ List *type_name_list = list_make2(makeString(nspname), makeString("scircle"));
+ TypeName *type_name = makeTypeNameFromNameList(type_name_list);
+ Oid type_oid = LookupTypeNameOid(NULL, type_name, true);
+ if (type_oid == InvalidOid)
+ elog(ERROR, "%s: unable to lookup type 'scircle'", __func__);
+ return type_oid;
+}
+
+static Oid
+scircleFunctionOid(Oid geotype, Oid callingfunc)
+{
+ const Oid radiustype = FLOAT8OID; /* Should always be FLOAT8OID */
+ const Oid scircle_function_args[2] = {geotype, radiustype};
+ const bool noError = true;
+ /* Expand function must be in same namespace as the caller */
+ char *nspname = get_namespace_name(get_func_namespace(callingfunc));
+ List *scircle_function_name = list_make2(makeString(nspname), makeString("scircle"));
+ Oid scircle_function_oid = LookupFuncName(scircle_function_name, 2, scircle_function_args, noError);
+ if (scircle_function_oid == InvalidOid)
+ elog(ERROR, "%s: unable to lookup 'scircle(Oid[%u], Oid[%u])'", __func__, geotype, radiustype);
+ return scircle_function_oid;
+}
+
+PG_FUNCTION_INFO_V1(spherepoint_dwithin_supportfn);
+Datum spherepoint_dwithin_supportfn(PG_FUNCTION_ARGS)
+{
+ Node *rawreq = (Node *) PG_GETARG_POINTER(0);
+ Node *ret = NULL;
+
+ if (IsA(rawreq, SupportRequestSelectivity))
+ {
+ SupportRequestSelectivity *req = (SupportRequestSelectivity *) rawreq;
+ Node *radiusarg = (Node *) list_nth(req->args, 2);
+ float8 selec;
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn sel called on %d", req->funcid)));
+
+ /*
+ * If the radius is a constant, compute the circle constant.
+ */
+ if (IsA(radiusarg, Const))
+ {
+ Const *constarg = (Const *) radiusarg;
+ float8 radius = DatumGetFloat8(constarg->constvalue);
+ selec = spherecircle_area_float(radius) / SPHERE_SURFACE;
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn const radius %g", radius)));
+ }
+ else
+ {
+ selec = DEFAULT_SCIRCLE_SEL;
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn non-const radius")));
+ }
+
+ if (req->is_join)
+ {
+ req->selectivity = selec;
+ }
+ else
+ {
+ req->selectivity = selec;
+ }
+ CLAMP_PROBABILITY(req->selectivity);
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn selectivity %g is_join %d", req->selectivity, req->is_join)));
+ ret = rawreq;
+ }
+
+ else if (IsA(rawreq, SupportRequestIndexCondition))
+ {
+ SupportRequestIndexCondition *req = (SupportRequestIndexCondition *) rawreq;
+
+ FuncExpr *clause = (FuncExpr *) req->node;
+ Oid funcid = clause->funcid;
+ Oid opfamilyoid = req->opfamily; /* OPERATOR FAMILY of the index */
+
+ Node *leftarg, *rightarg, *radiusarg;
+ Oid leftdatatype, oproid;
+
+ Oid scircle_type_oid = scircleTypeOid(funcid);
+ Expr *scircle_expr;
+ Expr *expr;
+
+ /*
+ * Extract "leftarg" as the arg matching the index and "rightarg" as
+ * the other, even if they were in the opposite order in the call.
+ */
+ if (req->indexarg == 0)
+ {
+ leftarg = linitial(clause->args);
+ rightarg = lsecond(clause->args);
+ }
+ else if (req->indexarg == 1)
+ {
+ rightarg = linitial(clause->args);
+ leftarg = lsecond(clause->args);
+ }
+ else
+ PG_RETURN_POINTER((Node *)NULL);
+
+ leftdatatype = exprType(leftarg);
+ Assert(leftdatatype == exprType(rightarg)); /* expect spoint, spoint */
+
+ radiusarg = (Node *) list_nth(clause->args, 2);
+
+ /*
+ * Given the index operator family and the arguments and the desired
+ * strategy number we can now lookup the operator we want.
+ */
+ oproid = get_opfamily_member(opfamilyoid,
+ leftdatatype,
+ scircle_type_oid,
+ 37); /* spoint <@ scircle */
+ if (!OidIsValid(oproid))
+ elog(ERROR,
+ "no spatial operator found for '%s': opfamily %u types %d %d strategy %d",
+ "scircle",
+ opfamilyoid,
+ leftdatatype,
+ scircle_type_oid,
+ 37);
+
+ /*
+ * If both the right argument and the radius are a constant, compute
+ * the circle constant. (makeFuncExpr won't constify by itself
+ * unfortunately.)
+ */
+ if (IsA(rightarg, Const) && IsA(radiusarg, Const))
+ {
+ Datum center = ((Const *) rightarg)->constvalue;
+ Datum radius = ((Const *) radiusarg)->constvalue;
+ Datum circle = DirectFunctionCall2(spherecircle_by_center, center, radius);
+ scircle_expr = (Expr *) makeConst(scircle_type_oid, -1, InvalidOid,
+ sizeof(SCIRCLE), circle, false, false);
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn index condition const")));
+ }
+ else
+ {
+ Oid scircle_function_oid = scircleFunctionOid(leftdatatype, clause->funcid);
+ scircle_expr = (Expr *) makeFuncExpr(scircle_function_oid, leftdatatype,
+ list_make2(rightarg, radiusarg),
+ InvalidOid, InvalidOid, COERCE_EXPLICIT_CALL);
+ ereport(DEBUG1, (errmsg("spherepoint_dwithin_supportfn index condition function")));
+ }
+
+ /*
+ * The comparison expression has to be a pseudo constant,
+ * (not volatile or dependent on the target index table)
+ */
+#if PG_VERSION_NUM >= 140000
+ if (!is_pseudo_constant_for_index(req->root, (Node*)scircle_expr, req->index))
+#else
+ if (!is_pseudo_constant_for_index((Node*)scircle_expr, req->index))
+#endif
+ PG_RETURN_POINTER((Node*)NULL);
+
+ /* OK, we can make an index expression */
+ expr = make_opclause(oproid, BOOLOID, false,
+ (Expr *) leftarg, scircle_expr,
+ InvalidOid, InvalidOid);
+
+ ret = (Node *)(list_make1(expr));
+
+ /* This is an exact index lookup */
+ req->lossy = false;
+ }
+
+ PG_RETURN_POINTER(ret);
+}
+
+#endif /* PG_VERSION_NUM */
diff --git a/src/pgs_util.h b/src/pgs_util.h
index b79170b..1b3ad76 100644
--- a/src/pgs_util.h
+++ b/src/pgs_util.h
@@ -14,6 +14,10 @@
#endif
#define EPSILON 1.0E-09 /* precision of floating point values */
+/* spherical circle constants */
+#define SPHERE_SURFACE (4 * PI)
+#define DEFAULT_SCIRCLE_SEL 1e-7
+
/* convert pg_sphere theta [pi/2 .. -pi/2] to healpix theta [0 .. pi] ([north .. south pole]) */
static inline double
conv_theta(double x)
@@ -31,4 +35,13 @@ static inline double deg_to_rad(double in)
return in * PI / 180;
}
+/*
+ * Area of circle on sphere
+ */
+static inline double
+spherecircle_area_float(double radius)
+{
+ return PID * (1 - cos(radius));
+}
+
#endif
diff --git a/src/point.c b/src/point.c
index 0cbfa00..498ee5e 100644
--- a/src/point.c
+++ b/src/point.c
@@ -8,6 +8,7 @@ PG_FUNCTION_INFO_V1(spherepoint_in);
PG_FUNCTION_INFO_V1(spherepoint_from_long_lat);
PG_FUNCTION_INFO_V1(spherepoint_from_long_lat_deg);
PG_FUNCTION_INFO_V1(spherepoint_distance);
+PG_FUNCTION_INFO_V1(spherepoint_dwithin);
PG_FUNCTION_INFO_V1(spherepoint_long);
PG_FUNCTION_INFO_V1(spherepoint_lat);
PG_FUNCTION_INFO_V1(spherepoint_x);
@@ -223,6 +224,17 @@ spherepoint_distance(PG_FUNCTION_ARGS)
}
+Datum
+spherepoint_dwithin(PG_FUNCTION_ARGS)
+{
+ SPoint *p1 = (SPoint *) PG_GETARG_POINTER(0);
+ SPoint *p2 = (SPoint *) PG_GETARG_POINTER(1);
+ float8 within = PG_GETARG_FLOAT8(2);
+ float8 dist = spoint_dist(p1, p2);
+
+ PG_RETURN_BOOL(FPle(dist, within));
+}
+
Datum
spherepoint_long(PG_FUNCTION_ARGS)
{