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
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
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
}
-LWGEOM *
+LWGEOM*
lwgeom_construct_empty(uint8_t type, int srid, char hasz, char hasm)
{
switch(type)
}
}
+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;
+ }
+}
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);
+}
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
$$ 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
#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);
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
#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; // <GeomCache>
+ GSERIALIZED* geom1; //
+ GSERIALIZED* geom2; //
+ size_t geom1_size; //
+ size_t geom2_size; //
+ int32 argnum; // </GeomCache>
+ CIRC_NODE* index;
+} CircTreeGeomCache;
+
+
+
/**
* Builder, freeer and public accessor for cached CIRC_NODE trees
*/
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;
+ }
+
+}
#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; // <GeomCache>
- GSERIALIZED* geom1; //
- GSERIALIZED* geom2; //
- size_t geom1_size; //
- size_t geom2_size; //
- int32 argnum; // </GeomCache>
- 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);
projects / postgis / commitdiff