From 0396a0468c349690aaed489794d4af829e310da0 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 1 Oct 2009 19:50:44 +0000 Subject: [PATCH] ST_Distance(geography, geography) roughed in. Small detail, currently returns answers in radians. :) git-svn-id: http://svn.osgeo.org/postgis/trunk@4573 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/lwgeodetic.c | 105 ++++++++++++++++++++++++++++++++++- postgis/Makefile.in | 3 +- postgis/geography.h | 31 +++++++++++ postgis/geography.sql.in.c | 55 +++++++++++------- postgis/geography_distance.c | 81 +++++++++++++++++++++++++++ postgis/geography_gist.c | 40 +++++-------- 6 files changed, 267 insertions(+), 48 deletions(-) create mode 100644 postgis/geography_distance.c diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 7ab69dc36..f4d38b885 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1348,14 +1348,117 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX *gbox1, GBO if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) || ( type2 == POLYGONTYPE && type1 == LINETYPE ) ) { + POINT2D p; + LWPOLY *lwpoly; + LWLINE *lwline; + GBOX *gbox; + double distance = MAXFLOAT; + int i; + + if( type1 == LINETYPE ) + { + lwline = (LWLINE*)lwgeom1; + lwpoly = (LWPOLY*)lwgeom2; + gbox = gbox2; + } + else + { + lwline = (LWLINE*)lwgeom2; + lwpoly = (LWPOLY*)lwgeom1; + gbox = gbox1; + } + getPoint2d_p(lwline->points, 0, &p); + + /* Point in polygon implies zero distance */ + if( lwpoly_covers_point2d(lwpoly, gbox, p) ) + return 0.0; + + /* Not contained, so what's the actual distance? */ + for( i = 0; i < lwpoly->nrings; i++ ) + { + double ring_distance = ptarray_distance_sphere(lwpoly->rings[i], lwline->points, tolerance); + if( ring_distance < distance ) + distance = ring_distance; + if( distance < tolerance ) + return distance; + } + return distance; + } /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */ if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) || ( type2 == POLYGONTYPE && type1 == LINETYPE ) ) { - } + POINT2D p; + LWPOLY *lwpoly1 = (LWPOLY*)lwgeom1; + LWPOLY *lwpoly2 = (LWPOLY*)lwgeom2; + double distance = MAXFLOAT; + int i, j; + /* Point of 2 in polygon 1 implies zero distance */ + getPoint2d_p(lwpoly1->rings[0], 0, &p); + if( lwpoly_covers_point2d(lwpoly2, gbox2, p) ) + return 0.0; + + /* Point of 1 in polygon 2 implies zero distance */ + getPoint2d_p(lwpoly2->rings[0], 0, &p); + if( lwpoly_covers_point2d(lwpoly1, gbox1, p) ) + return 0.0; + + /* Not contained, so what's the actual distance? */ + for( i = 0; i < lwpoly1->nrings; i++ ) + { + for( j = 0; j < lwpoly2->nrings; j++ ) + { + double ring_distance = ptarray_distance_sphere(lwpoly1->rings[i], lwpoly2->rings[j], tolerance); + if( ring_distance < distance ) + distance = ring_distance; + if( distance < tolerance ) + return distance; + } + } + return distance; + } + + /* Recurse into collections */ + if( lwgeom_contains_subgeoms(type1) ) + { + int i; + double distance = MAXFLOAT; + LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1; + + for( i = 0; i < col->ngeoms; i++ ) + { + double geom_distance = lwgeom_distance_sphere(col->geoms[i], lwgeom2, gbox1, gbox2, tolerance); + if( geom_distance < distance ) + distance = geom_distance; + if( distance < tolerance ) + return distance; + } + return distance; + } + + /* Recurse into collections */ + if( lwgeom_contains_subgeoms(type2) ) + { + int i; + double distance = MAXFLOAT; + LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2; + + for( i = 0; i < col->ngeoms; i++ ) + { + double geom_distance = lwgeom_distance_sphere(lwgeom1, col->geoms[i], gbox1, gbox2, tolerance); + if( geom_distance < distance ) + distance = geom_distance; + if( distance < tolerance ) + return distance; + } + return distance; + } + + + lwerror("arguments include unsupported geometry type (%s, %s)", lwgeom_typename(type1), lwgeom_typename(type1)); return -1.0; } diff --git a/postgis/Makefile.in b/postgis/Makefile.in index 3e2e1c9a1..69522c79e 100644 --- a/postgis/Makefile.in +++ b/postgis/Makefile.in @@ -51,7 +51,8 @@ PG_OBJS=lwgeom_pg.o \ lwgeom_rtree.o \ geography_inout.o \ geography_gist.o \ - geography_estimate.o + geography_estimate.o \ + geography_distance.o # Objects to build using PGXS OBJS=$(PG_OBJS) diff --git a/postgis/geography.h b/postgis/geography.h index ad177b387..8e1f3de39 100644 --- a/postgis/geography.h +++ b/postgis/geography.h @@ -9,7 +9,38 @@ * **********************************************************************/ + +/********************************************************************** +** GIDX structure. +** +** This is an n-dimensional (practically, the +** implementation is 2-4 dimensions) box used for index keys. The +** coordinates are anonymous, so we can have any number of dimensions. +** The sizeof a GIDX is 1 + 2 * ndims * sizeof(float). +** The order of values in a GIDX is +** xmin,xmax,ymin,ymax,zmin,zmax... +*/ +typedef struct +{ + int32 varsize; + float c[1]; +} GIDX; + +/* +** For some GiST support functions, it is more efficient to allocate +** memory for our GIDX off the stack and cast that memory into a GIDX. +** But, GIDX is variable length, what to do? We'll bake in the assumption +** that no GIDX is more than 4-dimensional for now, and ensure that much +** space is available. +** 4 bytes varsize + 4 dimensions * 2 ordinates * 4 bytes float size = 36 bytes +*/ +#define GIDX_MAX_SIZE 36 + + /* 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); +GIDX* gidx_new(int ndims); +void gbox_from_gidx(GIDX *gidx, GBOX *gbox); diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index a8b1ccc80..23acdbbcb 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -367,46 +367,63 @@ CREATE OR REPLACE FUNCTION ST_AsKML(int4, geography, int4) -- CREATE OR REPLACE FUNCTION _ST_AsGeoJson(int4, geography, int4, int4) - RETURNS text + RETURNS text AS 'MODULE_PATHNAME','geography_as_geojson' LANGUAGE 'C' IMMUTABLE STRICT; -- ST_AsGeoJson(geography, precision) / version=1 options=0 CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography, int4) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson(1, $1, $2, 0)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson(1, $1, $2, 0)' + LANGUAGE 'SQL' IMMUTABLE STRICT; -- ST_AsGeoJson(geography) / precision=15 version=1 options=0 CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson(1, $1, 15, 0)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson(1, $1, 15, 0)' + LANGUAGE 'SQL' IMMUTABLE STRICT; -- ST_AsGeoJson(version, geography) / precision=15 options=0 CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson($1, $2, 15, 0)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson($1, $2, 15, 0)' + LANGUAGE 'SQL' IMMUTABLE STRICT; -- ST_AsGeoJson(version, geography, precision) / options=0 CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson($1, $2, $3, 0)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson($1, $2, $3, 0)' + LANGUAGE 'SQL' IMMUTABLE STRICT; -- ST_AsGeoJson(geography, precision, options) / version=1 CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography, int4, int4) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson(1, $1, $2, $3)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson(1, $1, $2, $3)' + LANGUAGE 'SQL' IMMUTABLE STRICT; -- ST_AsGeoJson(version, geography, precision,options) CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4, int4) - RETURNS TEXT - AS 'SELECT _ST_AsGeoJson($1, $2, $3, $4)' - LANGUAGE 'SQL' IMMUTABLE STRICT; + RETURNS TEXT + AS 'SELECT _ST_AsGeoJson($1, $2, $3, $4)' + LANGUAGE 'SQL' IMMUTABLE STRICT; +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +-- Measurement Functions +-- Availability: 1.5.0 +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +-- Stop calculation and return once distance is less than tolerance +CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8) + RETURNS float8 + AS 'MODULE_PATHNAME','geography_distance_sphere' + LANGUAGE 'C' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION ST_Distance(geography, geography) + RETURNS float8 + AS 'SELECT _ST_Distance($1, $2, 0.0)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + + +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- COMMIT; diff --git a/postgis/geography_distance.c b/postgis/geography_distance.c new file mode 100644 index 000000000..81713818b --- /dev/null +++ b/postgis/geography_distance.c @@ -0,0 +1,81 @@ +/********************************************************************** + * $Id: geography_inout.c 4535 2009-09-28 18:16:21Z colivier $ + * + * PostGIS - Spatial Types for PostgreSQL + * Copyright 2009 Paul Ramsey + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include "postgres.h" + +#include "../postgis_config.h" + +#include +#include +#include +#include +#include + +#include "libgeom.h" /* For standard geometry types. */ +#include "lwgeom_pg.h" /* For debugging macros. */ +#include "geography.h" /* For utility functions. */ + +Datum geography_distance_sphere(PG_FUNCTION_ARGS); + + +/* +** geography_in(cstring) returns *GSERIALIZED +*/ +PG_FUNCTION_INFO_V1(geography_distance_sphere); +Datum geography_distance_sphere(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GBOX gbox1; + GBOX gbox2; + GIDX *gidx1 = gidx_new(3); + GIDX *gidx2 = gidx_new(3); + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + double tolerance; + double distance; + + /* We need the bounding boxes in case of polygon calculations, + which requires them to generate a stab-line to test point-in-polygon. */ + geography_datum_gidx(PG_GETARG_DATUM(0), gidx1); + geography_datum_gidx(PG_GETARG_DATUM(0), gidx2); + gbox_from_gidx(gidx1, &gbox1); + gbox_from_gidx(gidx2, &gbox2); + pfree(gidx1); + pfree(gidx2); + + /* Get our geometry objects loaded into memory. */ + g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + lwgeom1 = lwgeom_from_gserialized(g1); + lwgeom2 = lwgeom_from_gserialized(g2); + + /* Read our tolerance value. */ + tolerance = PG_GETARG_FLOAT8(2); + + /* Calculate the distance */ + distance = lwgeom_distance_sphere(lwgeom1, lwgeom2, &gbox1, &gbox2, tolerance); + + /* Something went wrong... should already be eloged */ + if( distance < 0.0 ) + { + elog(ERROR, "Error in geography_distance_sphere calculation."); + PG_RETURN_NULL(); + } + + /* Clean up, but not all the way to the point arrays */ + lwgeom_release(lwgeom1); + lwgeom_release(lwgeom2); + + PG_RETURN_FLOAT8(distance); + +} + diff --git a/postgis/geography_gist.c b/postgis/geography_gist.c index f7f12af1b..aa4225df0 100644 --- a/postgis/geography_gist.c +++ b/postgis/geography_gist.c @@ -33,33 +33,8 @@ #include "lwgeom_pg.h" /* For debugging macros. */ #include "geography.h" /* For utility functions. */ -/********************************************************************************* -** GIDX structure. -** -** This is an n-dimensional (practically, the -** implementation is 2-4 dimensions) box used for index keys. The -** coordinates are anonymous, so we can have any number of dimensions. -** The sizeof a GIDX is 1 + 2 * ndims * sizeof(float). -** The order of values in a GIDX is -** xmin,xmax,ymin,ymax,zmin,zmax... -*/ -typedef struct -{ - int32 varsize; - float c[1]; -} GIDX; -/* -** For some GiST support functions, it is more efficient to allocate -** memory for our GIDX off the stack and cast that memory into a GIDX. -** But, GIDX is variable length, what to do? We'll bake in the assumption -** that no GIDX is more than 4-dimensional for now, and ensure that much -** space is available. -** 4 bytes varsize + 4 dimensions * 2 ordinates * 4 bytes float size = 36 bytes -*/ -#define GIDX_MAX_SIZE 36 - /* ** When is a node split not so good? If more than 90% of the entries ** end up in one of the children. @@ -117,7 +92,7 @@ Datum geography_overlaps(PG_FUNCTION_ARGS); /* Allocates a new GIDX on the heap of the requested dimensionality */ -static GIDX* gidx_new(int ndims) +GIDX* gidx_new(int ndims) { size_t size = GIDX_SIZE(ndims); GIDX *g = (GIDX*)palloc(size); @@ -381,6 +356,17 @@ static GIDX* gidx_from_gbox(GBOX box) } +void gbox_from_gidx(GIDX *a, GBOX *gbox) +{ + gbox->xmin = (double)GIDX_GET_MIN(a,0); + gbox->ymin = (double)GIDX_GET_MIN(a,1); + gbox->zmin = (double)GIDX_GET_MIN(a,2); + gbox->xmax = (double)GIDX_GET_MAX(a,0); + gbox->ymax = (double)GIDX_GET_MAX(a,1); + gbox->zmax = (double)GIDX_GET_MAX(a,2); +} + + /* ** Overlapping GIDX box test. ** @@ -502,7 +488,7 @@ static bool gidx_equals(GIDX *a, GIDX *b) ** full geography and return the box based on that. If no box is available, ** return G_FAILURE, otherwise G_SUCCESS. */ -static int geography_datum_gidx(Datum geography_datum, GIDX *gidx) +int geography_datum_gidx(Datum geography_datum, GIDX *gidx) { GSERIALIZED *gpart; int result = G_SUCCESS; -- 2.50.1