From a6f362e30e22b74ca0589c16ac95518169b5b6fe Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 19 Jun 2012 03:42:30 +0000 Subject: [PATCH] Add in SQL binding for circ-tree cached searching. Currently in *Cached variants to allow comparisons. git-svn-id: http://svn.osgeo.org/postgis/trunk@9949 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/liblwgeom.h.in | 1 + liblwgeom/liblwgeom_internal.h | 7 ++ liblwgeom/lwcollection.c | 8 ++ liblwgeom/lwgeom.c | 31 ++++- liblwgeom/lwpoly.c | 8 +- liblwgeom/ptarray.c | 6 + postgis/geography.sql.in.c | 38 ++++++ postgis/geography_measurement.c | 138 ++++++++++++++++++++ postgis/geography_measurement_trees.c | 174 +++++++++++++++++++++++++- postgis/geography_measurement_trees.h | 22 +--- 10 files changed, 410 insertions(+), 23 deletions(-) diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index f8f5cfe89..325d40f04 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -1084,6 +1084,7 @@ extern void lwpoly_reverse(LWPOLY *poly); extern void lwtriangle_reverse(LWTRIANGLE *triangle); extern char* lwgeom_summary(const LWGEOM *lwgeom, int offset); extern char* lwpoint_to_latlon(const LWPOINT *p, const char *format); +extern int lwgeom_startpoint(const LWGEOM* lwgeom, POINT4D* pt); /** * Ensure the outer ring is clockwise oriented and all inner rings diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index 368bfcdd1..85c7f5bec 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -320,6 +320,13 @@ LWLINE *lwline_clone_deep(const LWLINE *lwgeom); LWPOLY *lwpoly_clone_deep(const LWPOLY *lwgeom); LWCOLLECTION *lwcollection_clone_deep(const LWCOLLECTION *lwgeom); +/* +* Startpoint +*/ +int lwpoly_startpoint(const LWPOLY* lwpoly, POINT4D* pt); +int ptarray_startpoint(const POINTARRAY* pa, POINT4D* pt); +int lwcollection_startpoint(const LWCOLLECTION* col, POINT4D* pt); + /* * Write into *ret the coordinates of the closest point on * segment A-B to the reference input point R diff --git a/liblwgeom/lwcollection.c b/liblwgeom/lwcollection.c index 2adeb0847..dedc02962 100644 --- a/liblwgeom/lwcollection.c +++ b/liblwgeom/lwcollection.c @@ -546,3 +546,11 @@ int lwcollection_allows_subtype(int collectiontype, int subtype) return LW_FALSE; } +int +lwcollection_startpoint(const LWCOLLECTION* col, POINT4D* pt) +{ + if ( col->ngeoms < 1 ) + return LW_FAILURE; + + return lwgeom_startpoint(col->geoms[0], pt); +} \ No newline at end of file diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c index dc5f47d83..60e055d05 100644 --- a/liblwgeom/lwgeom.c +++ b/liblwgeom/lwgeom.c @@ -1524,7 +1524,7 @@ lwgeom_affine(LWGEOM *geom, const AFFINE *affine) } -LWGEOM * +LWGEOM* lwgeom_construct_empty(uint8_t type, int srid, char hasz, char hasm) { switch(type) @@ -1554,4 +1554,33 @@ lwgeom_construct_empty(uint8_t type, int srid, char hasz, char hasm) } } +int +lwgeom_startpoint(const LWGEOM* lwgeom, POINT4D* pt) +{ + if ( ! lwgeom ) + return LW_FAILURE; + + switch( lwgeom->type ) + { + case POINTTYPE: + return ptarray_startpoint(((LWPOINT*)lwgeom)->point, pt); + case TRIANGLETYPE: + case CIRCSTRINGTYPE: + case LINETYPE: + return ptarray_startpoint(((LWLINE*)lwgeom)->points, pt); + case POLYGONTYPE: + return lwpoly_startpoint((LWPOLY*)lwgeom, pt); + case CURVEPOLYTYPE: + case COMPOUNDTYPE: + case MULTIPOINTTYPE: + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + return lwcollection_startpoint((LWCOLLECTION*)lwgeom, pt); + default: + lwerror("int: unsupported geometry type: %s", + lwtype_name(lwgeom->type)); + return LW_FAILURE; + } +} diff --git a/liblwgeom/lwpoly.c b/liblwgeom/lwpoly.c index 1ec3c7ff1..058974484 100644 --- a/liblwgeom/lwpoly.c +++ b/liblwgeom/lwpoly.c @@ -499,5 +499,11 @@ lwpoly_is_closed(const LWPOLY *poly) return LW_TRUE; } - +int +lwpoly_startpoint(const LWPOLY* poly, POINT4D* pt) +{ + if ( poly->nrings < 1 ) + return LW_FAILURE; + return ptarray_startpoint(poly->rings[0], pt); +} diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c index 4b9416c0b..eab44e4d3 100644 --- a/liblwgeom/ptarray.c +++ b/liblwgeom/ptarray.c @@ -1404,3 +1404,9 @@ ptarray_affine(POINTARRAY *pa, const AFFINE *a) LWDEBUG(3, "lwgeom_affine_ptarray end"); } + +int +ptarray_startpoint(const POINTARRAY* pa, POINT4D* pt) +{ + return getPoint4d_p(pa, 0, pt); +} \ No newline at end of file diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index f2c9de48c..d257bced2 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -541,6 +541,44 @@ CREATE OR REPLACE FUNCTION ST_DWithin(text, text, float8) $$ SELECT ST_DWithin($1::geometry, $2::geometry, $3); $$ LANGUAGE 'sql' IMMUTABLE; + +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +-- TEMPORARY TESTING FUNCTIONS FOR CACHED IMPLEMENTATIONS OF GEOGRAPHY +-- DISTANCE AND DWITHIN +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +CREATE OR REPLACE FUNCTION _ST_DistanceCached(geography, geography, float8, boolean) + RETURNS float8 + AS 'MODULE_PATHNAME','geography_distance_cached' + LANGUAGE 'c' IMMUTABLE STRICT + COST 100; + +CREATE OR REPLACE FUNCTION ST_DistanceCached(geography, geography, boolean) + RETURNS float8 + AS 'SELECT _ST_DistanceCached($1, $2, 0.0, $3)' + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION ST_DistanceCached(geography, geography) + RETURNS float8 + AS 'SELECT _ST_DistanceCached($1, $2, 0.0, true)' + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION _ST_DWithinCached(geography, geography, float8, boolean) + RETURNS boolean + AS 'MODULE_PATHNAME','geography_dwithin_cached' + LANGUAGE 'c' IMMUTABLE STRICT + COST 100; + +CREATE OR REPLACE FUNCTION ST_DWithinCached(geography, geography, float8) + RETURNS boolean + AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_DWithinCached($1, $2, $3, true)' + LANGUAGE 'sql' IMMUTABLE; +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +-- ---------- ---------- ---------- ---------- ---------- ---------- ---------- + + + + -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Area(geog geography, use_spheroid boolean DEFAULT true) RETURNS float8 diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 467c5e845..9d66a7fcd 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -24,10 +24,13 @@ #include "liblwgeom_internal.h" /* For FP comparators. */ #include "lwgeom_pg.h" /* For debugging macros. */ #include "geography.h" /* For utility functions. */ +#include "geography_measurement_trees.h" /* For circ_tree caching */ #include "lwgeom_transform.h" /* For SRID functions */ Datum geography_distance(PG_FUNCTION_ARGS); +Datum geography_distance_cached(PG_FUNCTION_ARGS); Datum geography_dwithin(PG_FUNCTION_ARGS); +Datum geography_dwithin_cached(PG_FUNCTION_ARGS); Datum geography_area(PG_FUNCTION_ARGS); Datum geography_length(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); @@ -99,6 +102,141 @@ Datum geography_distance(PG_FUNCTION_ARGS) PG_RETURN_FLOAT8(distance); } + +/* +** geography_distance_cached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid) +** returns double distance in meters +*/ +PG_FUNCTION_INFO_V1(geography_distance_cached); +Datum geography_distance_cached(PG_FUNCTION_ARGS) +{ + GSERIALIZED* g1 = NULL; + GSERIALIZED* g2 = NULL; + double distance; + double tolerance; + bool use_spheroid; + SPHEROID s; + + /* 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)); + + /* Read our tolerance value. */ + tolerance = PG_GETARG_FLOAT8(2); + + /* Read our calculation type. */ + use_spheroid = PG_GETARG_BOOL(3); + + /* Initialize spheroid */ + spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s); + + /* Set to sphere if requested */ + if ( ! use_spheroid ) + s.a = s.b = s.radius; + + /* Return NULL on empty arguments. */ + if ( gserialized_is_empty(g1) || gserialized_is_empty(g2) ) + { + PG_FREE_IF_COPY(g1, 0); + PG_FREE_IF_COPY(g2, 1); + PG_RETURN_NULL(); + } + + /* Do the brute force calculation if the cached calculation doesn't tick over */ + if ( LW_FAILURE == geography_distance_cache(fcinfo, g1, g2, &s, &distance) ) + { + LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1); + LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2); + distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, FP_TOLERANCE); + lwgeom_free(lwgeom1); + lwgeom_free(lwgeom2); + } + + /* Clean up */ + PG_FREE_IF_COPY(g1, 0); + PG_FREE_IF_COPY(g2, 1); + + /* Something went wrong, negative return... should already be eloged, return NULL */ + if ( distance < 0.0 ) + PG_RETURN_NULL(); + + PG_RETURN_FLOAT8(distance); +} + + +/* +** geography_dwithin(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid) +** returns double distance in meters +*/ +PG_FUNCTION_INFO_V1(geography_dwithin_cached); +Datum geography_dwithin_cached(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + double tolerance; + double distance; + bool use_spheroid; + SPHEROID s; + int dwithin = LW_FALSE; + + /* 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)); + + /* Read our tolerance value. */ + tolerance = PG_GETARG_FLOAT8(2); + + /* Read our calculation type. */ + use_spheroid = PG_GETARG_BOOL(3); + + /* Initialize spheroid */ + spheroid_init_from_srid(fcinfo, gserialized_get_srid(g1), &s); + + /* Set to sphere if requested */ + if ( ! use_spheroid ) + s.a = s.b = s.radius; + + /* Return FALSE on empty arguments. */ + if ( gserialized_is_empty(g1) || gserialized_is_empty(g2) ) + { + PG_FREE_IF_COPY(g1, 0); + PG_FREE_IF_COPY(g2, 1); + PG_RETURN_BOOL(FALSE); + } + + lwgeom1 = lwgeom_from_gserialized(g1); + lwgeom2 = lwgeom_from_gserialized(g2); + + /* Do the brute force calculation if the cached calculation doesn't tick over */ + if ( LW_FAILURE == geography_dwithin_cache(fcinfo, g1, g2, &s, tolerance, &dwithin) ) + { + LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1); + LWGEOM* lwgeom2 = lwgeom_from_gserialized(g2); + distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance); + dwithin = (distance <= tolerance); + lwgeom_free(lwgeom1); + lwgeom_free(lwgeom2); + } + + /* Clean up */ + lwgeom_free(lwgeom1); + lwgeom_free(lwgeom2); + PG_FREE_IF_COPY(g1, 0); + PG_FREE_IF_COPY(g2, 1); + + /* Something went wrong... should already be eloged, return FALSE */ + if ( distance < 0.0 ) + { + elog(ERROR, "lwgeom_distance_spheroid returned negative!"); + PG_RETURN_BOOL(FALSE); + } + + PG_RETURN_BOOL(dwithin); +} + + /* ** geography_dwithin(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid) ** returns double distance in meters diff --git a/postgis/geography_measurement_trees.c b/postgis/geography_measurement_trees.c index 9c9cc7fc5..e873f1547 100644 --- a/postgis/geography_measurement_trees.c +++ b/postgis/geography_measurement_trees.c @@ -1,5 +1,25 @@ #include "geography_measurement_trees.h" + +/* +* Specific tree types include all the generic slots and +* their own slots for their trees. We put the implementation +* for the CircTreeGeomCache here because we can't shove +* the PgSQL specific bits of the code (fcinfo) back into +* liblwgeom, where most of the circtree logic lives. +*/ +typedef struct { + int type; // + GSERIALIZED* geom1; // + GSERIALIZED* geom2; // + size_t geom1_size; // + size_t geom2_size; // + int32 argnum; // + CIRC_NODE* index; +} CircTreeGeomCache; + + + /** * Builder, freeer and public accessor for cached CIRC_NODE trees */ @@ -50,9 +70,161 @@ static GeomCacheMethods CircTreeCacheMethods = CircTreeAllocator }; -CircTreeGeomCache* +static CircTreeGeomCache* GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2) { return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2); } +static int +CircTreePIP(const CircTreeGeomCache* tree_cache, const GSERIALIZED* g, const LWGEOM* lwgeom) +{ + int tree_type = gserialized_get_type(g); + GBOX gbox; + GEOGRAPHIC_POINT gp; + POINT3D gp3; + POINT4D pt; + + if ( tree_type == POLYGONTYPE || tree_type == MULTIPOLYGONTYPE ) + { + /* Need a gbox to calculate an outside point */ + if ( LW_FAILURE == gserialized_get_gbox_p(g, &gbox) ) + { + LWGEOM* lwgeom_cached = lwgeom_from_gserialized(g); + lwgeom_calculate_gbox_geodetic(lwgeom_cached, &gbox); + lwgeom_free(lwgeom_cached); + } + + /* Need one point from the candidate geometry */ + if ( LW_FAILURE == lwgeom_startpoint(lwgeom, &pt) ) + { + lwerror("CircTreePIP unable to generate start point for lwgeom %p", lwgeom); + return LW_FALSE; + } + + /* Flip the candidate point into geographics */ + geographic_point_init(pt.x, pt.y, &gp); + geog2cart(&gp, &gp3); + + /* If the candidate isn't in the tree box, it's not in the tree area */ + if ( ! gbox_contains_point3d(&gbox, &gp3) ) + { + return LW_FALSE; + } + /* The candidate point is in the box, so it *might* be inside the tree */ + else + { + POINT2D pt_outside; /* latlon */ + POINT2D pt_inside; + pt_inside.x = pt.x; pt_inside.y = pt.y; + /* Calculate a definitive outside point */ + gbox_pt_outside(&gbox, &pt_outside); + /* Test the candidate point for strict containment */ + return circ_tree_contains_point(tree_cache->index, &pt_inside, &pt_outside, NULL); + } + + } + else + { + return LW_FALSE; + } +} + +int +geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance) +{ + CircTreeGeomCache* tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2); + + Assert(distance); + + if ( tree_cache && tree_cache->argnum && tree_cache->index ) + { + CIRC_NODE* circ_tree = NULL; + const GSERIALIZED* g = NULL; + LWGEOM* lwgeom = NULL; + + /* We need to dynamically build a tree for the uncached side of the function call */ + if ( tree_cache->argnum == 1 ) + { + lwgeom = lwgeom_from_gserialized(g2); + g = g1; + } + else if ( tree_cache->argnum == 2 ) + { + lwgeom = lwgeom_from_gserialized(g1); + g = g2; + } + else + lwerror("geography_distance_cache failed! This will never happen!"); + + if ( LW_TRUE == CircTreePIP(tree_cache, g, lwgeom) ) + { + *distance = 0.0; + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + + /* We do tree/tree distance, so turn the candidate geometry into a tree */ + circ_tree = lwgeom_calculate_circ_tree(lwgeom); + /* Calculate tree/tree distance */ + *distance = circ_tree_distance_tree(tree_cache->index, circ_tree, s, FP_TOLERANCE); + circ_tree_free(circ_tree); + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + else + { + return LW_FAILURE; + } +} + +int +geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin) +{ + CircTreeGeomCache* tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2); + double distance; + + Assert(dwithin); + + if ( tree_cache && tree_cache->argnum && tree_cache->index ) + { + CIRC_NODE* circ_tree = NULL; + const GSERIALIZED* g = NULL; + LWGEOM* lwgeom = NULL; + + /* We need to dynamically build a tree for the uncached side of the function call */ + if ( tree_cache->argnum == 1 ) + { + lwgeom = lwgeom_from_gserialized(g2); + g = g1; + } + else if ( tree_cache->argnum == 2 ) + { + lwgeom = lwgeom_from_gserialized(g1); + g = g2; + } + else + lwerror("geography_dwithin_cache failed! This will never happen!"); + + if ( LW_TRUE == CircTreePIP(tree_cache, g, lwgeom) ) + { + *dwithin = LW_TRUE; + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + + /* We do tree/tree distance, so turn the candidate geometry into a tree */ + circ_tree = lwgeom_calculate_circ_tree(lwgeom); + /* Calculate tree/tree distance */ + distance = circ_tree_distance_tree(tree_cache->index, circ_tree, s, tolerance); + *dwithin = (distance <= tolerance ? LW_TRUE : LW_FALSE); + circ_tree_free(circ_tree); + lwgeom_free(lwgeom); + return LW_SUCCESS; + } + else + { + return LW_FAILURE; + } + +} diff --git a/postgis/geography_measurement_trees.h b/postgis/geography_measurement_trees.h index 7ecd26183..3c7bc240b 100644 --- a/postgis/geography_measurement_trees.h +++ b/postgis/geography_measurement_trees.h @@ -2,23 +2,5 @@ #include "lwgeodetic_tree.h" #include "lwgeom_cache.h" -/* -* Specific tree types include all the generic slots and -* their own slots for their trees. We put the implementation -* for the CircTreeGeomCache here because we can't shove -* the PgSQL specific bits of the code (fcinfo) back into -* liblwgeom, where most of the circtree logic lives. -*/ -typedef struct { - int type; // - GSERIALIZED* geom1; // - GSERIALIZED* geom2; // - size_t geom1_size; // - size_t geom2_size; // - int32 argnum; // - CIRC_NODE* index; -} CircTreeGeomCache; - - - -CircTreeGeomCache* GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2); +int geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin); +int geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance); -- 2.40.0