Apply handling for EMPTY geometries to all geography functions per the DevWikiEmptyGe...

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 cafb3a1b6028ff5777947099a62b68fee2717d35..467975041fd2f8501015d00ed4bea9587665bd3e 100644 (file)
--- 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 cc0975e1d420105e6a2c3cdf416bc134937acc52..4f007dd84bb0de94170662ac812d349400166387 100644 (file)
--- 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 4c71287d00848fc26a9c3d9bd31f954d55cecb89..e4e9df28e90963c1931abe08b9d80fb135e3db9f 100644 (file)
--- 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 5eec1ad4750c8c3bcdef12dbcef840fcefb1df79..b16051867dbe85642c3523b8ea69ffc2ec9d1196 100644 (file)
--- 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 d8531be5c07a50992c6bc7f60e877957d8eeb281..3c54488bf5cd7b66f4b516ac27f75517f23df7c4 100644 (file)
--- 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);