From: Paul Ramsey Date: Tue, 10 Nov 2009 19:30:43 +0000 (+0000) Subject: Apply handling for EMPTY geometries to all geography functions per the DevWikiEmptyGe... X-Git-Tag: 1.5.0b1~263 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=2d305b179574388a4ecb3bd127f22f3ee02950ed;p=postgis Apply handling for EMPTY geometries to all geography functions per the DevWikiEmptyGeometry page. git-svn-id: http://svn.osgeo.org/postgis/trunk@4778 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index cafb3a1b6..467975041 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1742,6 +1742,7 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) * Calculate the distance between two LWGEOMs, using the coordinates are * longitude and latitude. Return immediately when the calulated distance drops * below the tolerance (useful for dwithin calculations). +* Return a negative distance for incalculable cases. */ double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, double tolerance) { @@ -1753,10 +1754,11 @@ double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GB LWDEBUGF(4, "entered function, tolerance %.8g", tolerance); - /* What's the distance to an empty geometry? We don't know. */ + /* What's the distance to an empty geometry? We don't know. + Return a negative number so the caller can catch this case. */ if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) { - return 0.0; + return -1.0; } type1 = TYPE_GETTYPE(lwgeom1->type); diff --git a/postgis/geography.h b/postgis/geography.h index cc0975e1d..4f007dd84 100644 --- a/postgis/geography.h +++ b/postgis/geography.h @@ -86,6 +86,8 @@ typedef struct 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); +/* Pull out the gidx bounding box from an already de-toasted geography */ +int geography_gidx(GSERIALIZED *g, GIDX *gidx); /* Convert a gidx to a gbox */ void gbox_from_gidx(GIDX *gidx, GBOX *gbox); /* Convert a gbox to a new gidx */ diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 4c71287d0..e4e9df28e 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -471,6 +471,14 @@ CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8, boolean) LANGUAGE 'C' IMMUTABLE STRICT COST 100; +-- Stop calculation and return immediately once distance is less than tolerance +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION _ST_DWithin(geography, geography, float8, boolean) + RETURNS boolean + AS 'MODULE_PATHNAME','geography_dwithin' + LANGUAGE 'C' IMMUTABLE STRICT + COST 100; + -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Distance(geography, geography, boolean) RETURNS float8 @@ -501,14 +509,14 @@ CREATE OR REPLACE FUNCTION _ST_Expand(geography, float8) -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8, boolean) RETURNS boolean - AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3, $4) < $3' + AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_DWithin($1, $2, $3, $4)' LANGUAGE 'SQL' IMMUTABLE; -- Currently defaulting to spheroid calculations -- 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, true) < $3' + AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_DWithin($1, $2, $3, true)' LANGUAGE 'SQL' IMMUTABLE; -- Availability: 1.5.0 - this is just a hack to prevent unknown from causing ambiguous name because of geography diff --git a/postgis/geography_gist.c b/postgis/geography_gist.c index 5eec1ad47..b16051867 100644 --- a/postgis/geography_gist.c +++ b/postgis/geography_gist.c @@ -586,6 +586,46 @@ int geography_datum_gidx(Datum geography_datum, GIDX *gidx) return result; } +/* +** Peak into a geography to find the bounding box. If the +** box is there, copy it out and return it. If not, calculate the box from the +** full geography and return the box based on that. If no box is available, +** return G_FAILURE, otherwise G_SUCCESS. +*/ +int geography_gidx(GSERIALIZED *g, GIDX *gidx) +{ + int result = G_SUCCESS; + + POSTGIS_DEBUG(4, "entered function"); + + POSTGIS_DEBUGF(4, "got flags %d", gpart->flags); + + if ( FLAGS_GET_BBOX(g->flags) && FLAGS_GET_GEODETIC(g->flags) ) + { + const size_t size = 2 * 3 * sizeof(float); + POSTGIS_DEBUG(4, "copying box out of serialization"); + memcpy(gidx->c, g->data, size); + SET_VARSIZE(gidx, VARHDRSZ + size); + } + else + { + GBOX gbox; + POSTGIS_DEBUG(4, "calculating new box from scratch"); + if( gserialized_calculate_gbox_geocentric_p(g, &gbox) == G_FAILURE ) + { + POSTGIS_DEBUG(4, "calculated null bbox, returning null"); + return G_FAILURE; + } + result = gidx_from_gbox_p(gbox, gidx); + } + if( result == G_SUCCESS ) + { + POSTGIS_DEBUGF(4, "got gidx %s", gidx_to_string(gidx)); + } + + return result; +} + /*********************************************************************** * GiST Support Functions */ @@ -603,11 +643,14 @@ Datum geography_overlaps(PG_FUNCTION_ARGS) GIDX *gbox1 = (GIDX*)gboxmem1; GIDX *gbox2 = (GIDX*)gboxmem2; - geography_datum_gidx(PG_GETARG_DATUM(0), gbox1); - geography_datum_gidx(PG_GETARG_DATUM(1), gbox2); - - if ( gidx_overlaps(gbox1, gbox2) == LW_TRUE ) + /* Must be able to build box for each arguement (ie, not empty geometry) + and overlap boxes to return true. */ + if( geography_datum_gidx(PG_GETARG_DATUM(0), gbox1) && + geography_datum_gidx(PG_GETARG_DATUM(1), gbox2) && + gidx_overlaps(gbox1, gbox2) ) + { PG_RETURN_BOOL(TRUE); + } PG_RETURN_BOOL(FALSE); } diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index d8531be5c..3c54488bf 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -24,6 +24,7 @@ #include "geography.h" /* For utility functions. */ Datum geography_distance(PG_FUNCTION_ARGS); +Datum geography_dwithin(PG_FUNCTION_ARGS); Datum geography_area(PG_FUNCTION_ARGS); Datum geography_length(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); @@ -32,7 +33,7 @@ Datum geography_covers(PG_FUNCTION_ARGS); Datum geography_bestsrid(PG_FUNCTION_ARGS); /* -** geography_distance_sphere(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance) +** geography_distance(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid) ** returns double distance in meters */ PG_FUNCTION_INFO_V1(geography_distance); @@ -44,18 +45,18 @@ Datum geography_distance(PG_FUNCTION_ARGS) GBOX gbox2; GSERIALIZED *g1 = NULL; GSERIALIZED *g2 = NULL; - double tolerance; 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); @@ -66,32 +67,107 @@ Datum geography_distance(PG_FUNCTION_ARGS) if( ! use_spheroid ) s.a = s.b = s.radius; + lwgeom1 = lwgeom_from_gserialized(g1); + lwgeom2 = lwgeom_from_gserialized(g2); + + /* Return NULL on empty arguments. */ + if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) + { + PG_RETURN_NULL(); + } + /* We need the bounding boxes in case of polygon calculations, which requires them to generate a stab-line to test point-in-polygon. */ if( ! gbox_from_gserialized(g1, &gbox1) || ! gbox_from_gserialized(g2, &gbox2) ) { - elog(ERROR, "Error in gbox_from_gserialized calculation."); + elog(NOTICE, "gbox_from_gserialized unable to calculate bounding box!"); + PG_RETURN_NULL(); + } + + distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, s, 0.0); + + /* Something went wrong, negative return... should already be eloged, return NULL */ + if( distance < 0.0 ) + { PG_RETURN_NULL(); } + + /* Clean up, but not all the way to the point arrays */ + lwgeom_release(lwgeom1); + lwgeom_release(lwgeom2); + + 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); +Datum geography_dwithin(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GBOX gbox1; + GBOX gbox2; + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + double tolerance; + double distance; + 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(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + /* Set to sphere if requested */ + if( ! use_spheroid ) + s.a = s.b = s.radius; + lwgeom1 = lwgeom_from_gserialized(g1); lwgeom2 = lwgeom_from_gserialized(g2); + /* Return FALSE on empty arguments. */ + if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) + { + PG_RETURN_BOOL(FALSE); + } + + /* We need the bounding boxes in case of polygon calculations, + which requires them to generate a stab-line to test point-in-polygon. */ + if( ! gbox_from_gserialized(g1, &gbox1) || + ! gbox_from_gserialized(g2, &gbox2) ) + { + elog(NOTICE, "gbox_from_gserialized unable to calculate bounding box!"); + PG_RETURN_BOOL(FALSE); + } + distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, s, tolerance); - /* Something went wrong... should already be eloged */ + /* Something went wrong... should already be eloged, return FALSE */ if( distance < 0.0 ) { - elog(ERROR, "lwgeom_distance_spheroid returned < 0.0"); - PG_RETURN_NULL(); + elog(ERROR, "lwgeom_distance_spheroid returned negative!"); + PG_RETURN_BOOL(FALSE); } /* Clean up, but not all the way to the point arrays */ lwgeom_release(lwgeom1); lwgeom_release(lwgeom2); - PG_RETURN_FLOAT8(distance); + PG_RETURN_BOOL(distance < tolerance); } @@ -113,12 +189,19 @@ Datum geography_expand(PG_FUNCTION_ARGS) 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)); + + /* Get our bounding box out of the geography, return right away if + input is an EMPTY geometry. */ + if ( geography_gidx(g, gidx) == G_FAILURE ) + { + g_out = palloc(VARSIZE(g)); + memcpy(g_out, g, VARSIZE(g)); + pfree(gidx); + PG_RETURN_POINTER(g_out); + } /* Read our distance value and normalize to unit-sphere. */ distance = PG_GETARG_FLOAT8(1) / WGS84_RADIUS; @@ -134,10 +217,12 @@ Datum geography_expand(PG_FUNCTION_ARGS) pfree(gidx); if( g_out == NULL ) + { + elog(ERROR, "gidx_insert_into_gserialized tried to insert mismatched dimensionality box into geography"); PG_RETURN_NULL(); + } PG_RETURN_POINTER(g_out); - } /* @@ -166,7 +251,16 @@ Datum geography_area(PG_FUNCTION_ARGS) /* User requests spherical calculation, turn our spheroid into a sphere */ if( ! use_spheroid ) s.a = s.b = s.radius; + + lwgeom = lwgeom_from_gserialized(g); + /* EMPTY things have no area */ + if( lwgeom_is_empty(lwgeom) ) + { + lwgeom_release(lwgeom); + PG_RETURN_FLOAT8(0.0); + } + /* We need the bounding box to get an outside point for area algorithm */ if( ! gbox_from_gserialized(g, &gbox) ) { @@ -174,8 +268,6 @@ Datum geography_area(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } - lwgeom = lwgeom_from_gserialized(g); - /* Test for cases that are currently not handled by spheroid code */ if ( use_spheroid ) { @@ -225,6 +317,13 @@ Datum geography_length(PG_FUNCTION_ARGS) g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); lwgeom = lwgeom_from_gserialized(g); + /* EMPTY things have no length */ + if( lwgeom_is_empty(lwgeom) ) + { + lwgeom_release(lwgeom); + PG_RETURN_FLOAT8(0.0); + } + /* Read our calculation type */ use_spheroid = PG_GETARG_BOOL(1); @@ -320,6 +419,18 @@ Datum geography_covers(PG_FUNCTION_ARGS) elog(ERROR, "geography_covers: only POLYGON and POINT types are currently supported"); PG_RETURN_NULL(); } + + /* Construct our working geometries */ + lwgeom1 = lwgeom_from_gserialized(g1); + lwgeom2 = lwgeom_from_gserialized(g2); + + /* EMPTY never intersects with another geometry */ + if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) + { + lwgeom_release(lwgeom1); + lwgeom_release(lwgeom2); + PG_RETURN_BOOL(false); + } /* We need the bounding boxes in case of polygon calculations, which requires them to generate a stab-line to test point-in-polygon. */ @@ -330,10 +441,6 @@ Datum geography_covers(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } - /* Construct our working geometries */ - lwgeom1 = lwgeom_from_gserialized(g1); - lwgeom2 = lwgeom_from_gserialized(g2); - /* Calculate answer */ result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2, gbox1, gbox2);