From: Paul Ramsey Date: Fri, 2 Oct 2009 18:28:06 +0000 (+0000) Subject: Largely untested implementation of ST_DWithin(geography, geography). X-Git-Tag: 1.5.0b1~431 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=3f9570136ed837695056fd523b56724a87aa8e09;p=postgis Largely untested implementation of ST_DWithin(geography, geography). We're in business baby! git-svn-id: http://svn.osgeo.org/postgis/trunk@4579 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/postgis/geography.h b/postgis/geography.h index 4aae1c6eb..931431be5 100644 --- a/postgis/geography.h +++ b/postgis/geography.h @@ -9,13 +9,29 @@ * **********************************************************************/ -/* + +/********************************************************************** ** Spherical radius. ** Moritz, H. (1980). Geodetic Reference System 1980, by resolution of the XVII General Assembly of the IUGG in Canberra. ** http://en.wikipedia.org/wiki/Earth_radius */ #define WGS84_RADIUS 6371009.0 +/* For reference, the old value used in distance_sphere was 6370986.884258304 */ + + +/********************************************************************** +** Useful functions for all GEOGRAPHY handlers. +*/ + +/* Convert lwgeom to newly allocated gserialized */ +GSERIALIZED* geography_serialize(LWGEOM *lwgeom); +/* Check that the typmod matches the flags on the lwgeom */ +void geography_valid_typmod(LWGEOM *lwgeom, int32 typmod); +/* Check that the type is legal in geography (no curves please!) */ +void geography_valid_type(uchar type); + + /********************************************************************** ** GIDX structure. @@ -43,11 +59,37 @@ typedef struct */ #define GIDX_MAX_SIZE 36 +/********************************************************************************* +** GIDX support functions. +** +** We put the GIDX support here rather than libgeom because it is a specialized +** type only used for indexing purposes. It also makes use of some PgSQL +** infrastructure like the VARSIZE() macros. +*/ -/* Useful functions for all GEOGRAPHY handlers. */ -GSERIALIZED* geography_serialize(LWGEOM *lwgeom); -void geography_valid_typmod(LWGEOM *lwgeom, int32 typmod); -void geography_valid_type(uchar type); -int geography_datum_gidx(Datum geography_datum, GIDX *gidx); +/* Returns number of dimensions for this GIDX */ +#define GIDX_NDIMS(gidx) ((VARSIZE((gidx)) - VARHDRSZ) / (2 * sizeof(float))) +/* Minimum accessor. */ +#define GIDX_GET_MIN(gidx, dimension) ((gidx)->c[2*(dimension)]) +/* Maximum accessor. */ +#define GIDX_GET_MAX(gidx, dimension) ((gidx)->c[2*(dimension)+1]) +/* Minimum setter. */ +#define GIDX_SET_MIN(gidx, dimension, value) ((gidx)->c[2*(dimension)] = (value)) +/* Maximum setter. */ +#define GIDX_SET_MAX(gidx, dimension, value) ((gidx)->c[2*(dimension)+1] = (value)) +/* Returns the size required to store a GIDX of requested dimension */ +#define GIDX_SIZE(dimensions) (sizeof(int32) + 2*(dimensions)*sizeof(float)) +/* Allocate a new gidx */ GIDX* gidx_new(int ndims); +/* Pull out the gidx bounding box with a absolute minimum system overhead */ +int geography_datum_gidx(Datum geography_datum, GIDX *gidx); +/* Convert a gidx to a gbox */ void gbox_from_gidx(GIDX *gidx, GBOX *gbox); +/* Convert a gbox to a new gidx */ +GIDX* gidx_from_gbox(GBOX box); +/* Copy a new bounding box into an existing gserialized */ +GSERIALIZED* gidx_insert_into_gserialized(GSERIALIZED *g, GIDX *gidx); + + + + diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 23acdbbcb..774d61549 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -412,18 +412,32 @@ CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4, int4) -- Availability: 1.5.0 -- ---------- ---------- ---------- ---------- ---------- ---------- ---------- --- Stop calculation and return once distance is less than tolerance +-- Stop calculation and return immediately once distance is less than tolerance +-- Availability: 1.5.0 CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8) RETURNS float8 AS 'MODULE_PATHNAME','geography_distance_sphere' LANGUAGE 'C' IMMUTABLE STRICT; +-- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Distance(geography, geography) RETURNS float8 AS 'SELECT _ST_Distance($1, $2, 0.0)' LANGUAGE 'SQL' IMMUTABLE STRICT; +-- Only expands the bounding box, the actual geometry will remain unchanged, use with care. +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION _ST_Expand(geography, float8) + RETURNS geography + AS 'MODULE_PATHNAME','geography_expand' + LANGUAGE 'C' IMMUTABLE STRICT; +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8) + RETURNS boolean + AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3) < $3' + LANGUAGE 'SQL' IMMUTABLE; + -- ---------- ---------- ---------- ---------- ---------- ---------- ---------- COMMIT; diff --git a/postgis/geography_distance.c b/postgis/geography_distance.c index 7a6e4b538..5a81cecfd 100644 --- a/postgis/geography_distance.c +++ b/postgis/geography_distance.c @@ -24,10 +24,12 @@ #include "geography.h" /* For utility functions. */ Datum geography_distance_sphere(PG_FUNCTION_ARGS); +Datum geography_expand(PG_FUNCTION_ARGS); /* -** geography_in(cstring) returns *GSERIALIZED +** geography_distance_sphere(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance) +** returns double distance in meters */ PG_FUNCTION_INFO_V1(geography_distance_sphere); Datum geography_distance_sphere(PG_FUNCTION_ARGS) @@ -72,6 +74,8 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } + /* Currently normalizing with a fixed WGS84 radius, in future this + should be the average radius of the SRID in play */ distance = distance * WGS84_RADIUS; /* Clean up, but not all the way to the point arrays */ @@ -82,3 +86,48 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) } + +/* +** geography_expand(GSERIALIZED *g) returns *GSERIALIZED +** +** warning, this tricky little function does not expand the +** geometry at all, just re-writes bounding box value to be +** a bit bigger. only useful when passing the result along to +** an index operator (&&) +*/ +PG_FUNCTION_INFO_V1(geography_expand); +Datum geography_expand(PG_FUNCTION_ARGS) +{ + GIDX *gidx = gidx_new(3); + GSERIALIZED *g = NULL; + GSERIALIZED *g_out = NULL; + double distance; + float fdistance; + int i; + + /* Get our bounding box out of the geography */ + geography_datum_gidx(PG_GETARG_DATUM(0), gidx); + + /* Get a pointer to the geography */ + g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + /* Read our distance value and normalize to unit-sphere. */ + distance = PG_GETARG_FLOAT8(1) / WGS84_RADIUS; + fdistance = (float)distance; + + for( i = 0; i < 3; i++ ) + { + GIDX_SET_MIN(gidx, i, GIDX_GET_MIN(gidx, i) - fdistance); + GIDX_SET_MAX(gidx, i, GIDX_GET_MAX(gidx, i) + fdistance); + } + + g_out = gidx_insert_into_gserialized(g, gidx); + pfree(gidx); + + if( g_out == NULL ) + PG_RETURN_NULL(); + + PG_RETURN_POINTER(g_out); + +} + diff --git a/postgis/geography_gist.c b/postgis/geography_gist.c index aa4225df0..c11e72dc9 100644 --- a/postgis/geography_gist.c +++ b/postgis/geography_gist.c @@ -344,7 +344,7 @@ static int gidx_from_gbox_p(GBOX box, GIDX *a) return G_SUCCESS; } -static GIDX* gidx_from_gbox(GBOX box) +GIDX* gidx_from_gbox(GBOX box) { int ndims; GIDX *a; @@ -366,6 +366,54 @@ void gbox_from_gidx(GIDX *a, GBOX *gbox) gbox->zmax = (double)GIDX_GET_MAX(a,2); } +/* +** Make a copy of a GSERIALIZED, with a new bounding box value embedded. +*/ +GSERIALIZED* gidx_insert_into_gserialized(GSERIALIZED *g, GIDX *gidx) +{ + int g_ndims = (FLAGS_GET_GEODETIC(g->flags) ? 3 : FLAGS_NDIMS(g->flags)); + int box_ndims = GIDX_NDIMS(gidx); + GSERIALIZED *g_out = NULL; + size_t box_size = 2 * g_ndims * sizeof(float); + + /* The dimensionality of the inputs has to match or we are SOL. */ + if( g_ndims != box_ndims ) + { + return NULL; + } + + /* Serialized already has room for a box. We just need to copy it and + write the new values into place. */ + if ( FLAGS_GET_BBOX(g->flags) ) + { + g_out = palloc(VARSIZE(g)); + memcpy(g_out, g, VARSIZE(g)); + } + /* Serialized has no box. We need to allocate enough space for the old + data plus the box, and leave a gap in the memory segment to write + the new values into. + */ + else + { + size_t varsize_new = VARSIZE(g) + box_size; + uchar *ptr; + g_out = palloc(varsize_new); + /* Copy the head of g into place */ + memcpy(g_out, g, 8); + /* Copy the body of g into place after leaving space for the box */ + ptr = g_out->data; + ptr += box_size; + memcpy(ptr, g->data, VARSIZE(g) - 8); + g_out->flags = FLAGS_SET_BBOX(g_out->flags, 1); + SET_VARSIZE(g_out, varsize_new); + } + + /* Now write the gidx values into the memory segement */ + memcpy(g_out->data, gidx->c, box_size); + + return g_out; +} + /* ** Overlapping GIDX box test. diff --git a/postgis/geography_inout.c b/postgis/geography_inout.c index 1840a1782..de28a26b5 100644 --- a/postgis/geography_inout.c +++ b/postgis/geography_inout.c @@ -438,7 +438,7 @@ Datum geography_typmod_out(PG_FUNCTION_ARGS) } /* -** geography_as_text(*GSERIALIZED) returns cstring +** geography_as_text(*GSERIALIZED) returns text */ PG_FUNCTION_INFO_V1(geography_as_text); Datum geography_as_text(PG_FUNCTION_ARGS) @@ -485,7 +485,7 @@ Datum geography_as_text(PG_FUNCTION_ARGS) /* -** geography_as_gml(*GSERIALIZED) returns cstring +** geography_as_gml(*GSERIALIZED) returns text */ PG_FUNCTION_INFO_V1(geography_as_gml); Datum geography_as_gml(PG_FUNCTION_ARGS) @@ -558,7 +558,7 @@ Datum geography_as_gml(PG_FUNCTION_ARGS) /* -** geography_as_kml(*GSERIALIZED) returns cstring +** geography_as_kml(*GSERIALIZED) returns text */ PG_FUNCTION_INFO_V1(geography_as_kml); Datum geography_as_kml(PG_FUNCTION_ARGS) @@ -614,7 +614,7 @@ Datum geography_as_kml(PG_FUNCTION_ARGS) /* -** geography_as_svg(*GSERIALIZED) returns cstring +** geography_as_svg(*GSERIALIZED) returns text */ PG_FUNCTION_INFO_V1(geography_as_svg); Datum geography_as_svg(PG_FUNCTION_ARGS) @@ -661,7 +661,7 @@ Datum geography_as_svg(PG_FUNCTION_ARGS) /* -** geography_as_geojson(*GSERIALIZED) returns cstring +** geography_as_geojson(*GSERIALIZED) returns text */ PG_FUNCTION_INFO_V1(geography_as_geojson); Datum geography_as_geojson(PG_FUNCTION_ARGS) @@ -767,7 +767,7 @@ Datum geography_from_text(PG_FUNCTION_ARGS) } /* -** geography_as_binary(*GSERIALIZED) returns *char +** geography_as_binary(*GSERIALIZED) returns bytea */ PG_FUNCTION_INFO_V1(geography_as_binary); Datum geography_as_binary(PG_FUNCTION_ARGS)