Fixed behavior of ST_ConvexHull() for empty rasters. Ticket #2126

git-svn-id: http://svn.osgeo.org/postgis/trunk@10810 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/NEWS b/NEWS
index 6652b19c647ab8a0ce3daf527aa2dd4525694f62..1482a756f4ea3c6b3c739a8f1c84449c5895d927 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -121,6 +121,7 @@ PostGIS 2.1.0
   - #2077, Fixed incorrect values returning from ST_Hillshade()
   - #2019, ST_FlipCoordinates does not update bbox
   - #2100, ST_AsRaster may not return raster with specified pixel type
+  - #2126, Better handling of empty rasters from ST_ConvexHull()
 
 PostGIS 2.0.3
 2013/MM/DD

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index 6ef2f37746718268dc6c31e15017f148223ea501..06a3a2a3ff31e29cf889e703d97d2a03886103a5 100644 (file)
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -6281,69 +6281,132 @@ rt_raster_gdal_polygonize(
        return pols;
 }
 
-LWPOLY*
-rt_raster_get_convex_hull(rt_raster raster) {
-               double gt[6] = {0.0};
-    POINTARRAY **rings = NULL;
-    POINTARRAY *pts = NULL;
-    LWPOLY* ret = NULL;
-    POINT4D p4d;
+/**
+ * Get raster's convex hull.
+ *
+ * The convex hull is typically a 4 vertices (5 to be closed)
+ * single ring polygon bearing the raster's rotation and using
+ * projection coordinates.
+ *
+ * @param raster : the raster to get info from
+ * @param **hull : pointer to convex hull
+ *
+ * @return ES_NONE if success, ES_ERROR if error
+ */
+rt_errorstate
+rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull) {
+       double gt[6] = {0.0};
+       int srid = SRID_UNKNOWN;
 
-    assert(NULL != raster);
+       POINTARRAY *pts = NULL;
+       POINT4D p4d;
 
-    RASTER_DEBUGF(3, "rt_raster_get_convex_hull: raster is %dx%d",
-            raster->width, raster->height);
+       assert(hull != NULL);
+       *hull = NULL;
 
-    if ((!raster->width) || (!raster->height)) {
-        return 0;
-    }
+       /* raster is NULL, convex hull is NULL */
+       if (raster == NULL)
+               return ES_NONE;
 
-    rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
-    if (!rings) {
-        rterror("rt_raster_get_convex_hull: Out of memory [%s:%d]", __FILE__, __LINE__);
-        return 0;
-    }
-    rings[0] = ptarray_construct(0, 0, 5);
-    /* TODO: handle error on ptarray construction */
-    /* XXX jorgearevalo: the error conditions aren't managed in ptarray_construct */
-    if (!rings[0]) {
-        rterror("rt_raster_get_convex_hull: Out of memory [%s:%d]", __FILE__, __LINE__);
-        return 0;
-    }
-    pts = rings[0];
-
-    /* Upper-left corner (first and last points) */
-    rt_raster_cell_to_geopoint(raster,
-            0, 0,
-            &p4d.x, &p4d.y,
-                                               gt);
-    ptarray_set_point4d(pts, 0, &p4d);
-    ptarray_set_point4d(pts, 4, &p4d); /* needed for closing it? */
-
-    /* Upper-right corner (we go clockwise) */
-    rt_raster_cell_to_geopoint(raster,
-            raster->width, 0,
-            &p4d.x, &p4d.y,
-                                               gt);
-    ptarray_set_point4d(pts, 1, &p4d);
-
-    /* Lower-right corner */
-    rt_raster_cell_to_geopoint(raster,
-            raster->width, raster->height,
-            &p4d.x, &p4d.y,
-                                               gt);
-    ptarray_set_point4d(pts, 2, &p4d);
-
-    /* Lower-left corner */
-    rt_raster_cell_to_geopoint(raster,
-            0, raster->height,
-            &p4d.x, &p4d.y,
-                                               gt);
-    ptarray_set_point4d(pts, 3, &p4d);
-
-    ret = lwpoly_construct(rt_raster_get_srid(raster), 0, 1, rings);
+       /* raster metadata */
+       srid = rt_raster_get_srid(raster);
+       rt_raster_get_geotransform_matrix(raster, gt);
 
-    return ret;
+       RASTER_DEBUGF(3, "rt_raster_get_convex_hull: raster is %dx%d", raster->width, raster->height); 
+
+       /* return point or line since at least one of the two dimensions is 0 */
+       if ((!raster->width) || (!raster->height)) {
+               p4d.x = gt[0];
+               p4d.y = gt[3];
+
+               /* return point */
+               if (!raster->width && !raster->height) {
+                       LWPOINT *point = lwpoint_make2d(srid, p4d.x, p4d.y);
+                       *hull = lwpoint_as_lwgeom(point);
+               }
+               /* return linestring */
+               else {
+                       LWLINE *line = NULL;
+                       pts = ptarray_construct_empty(0, 0, 2);
+
+                       /* first point of line */
+                       ptarray_append_point(pts, &p4d, LW_TRUE);
+
+                       /* second point of line */
+                       if (rt_raster_cell_to_geopoint(
+                               raster,
+                               rt_raster_get_width(raster), rt_raster_get_height(raster),
+                               &p4d.x, &p4d.y,
+                               gt
+                       ) != ES_NONE) {
+                               rterror("rt_raster_get_convex_hull: Unable to get second point for linestring");
+                               return ES_ERROR;
+                       }
+                       ptarray_append_point(pts, &p4d, LW_TRUE);
+                       line = lwline_construct(srid, NULL, pts);
+
+                       *hull = lwline_as_lwgeom(line);
+               }
+
+               return ES_NONE;
+       }
+       else {
+               POINTARRAY **rings = NULL;
+               LWPOLY* poly = NULL;
+
+               /* only one ring */
+               rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
+               if (!rings) {
+                       rterror("rt_raster_get_convex_hull: Unable to allocate memory for polygon ring");
+                       return ES_ERROR;
+               }
+               rings[0] = ptarray_construct(0, 0, 5);
+               /* TODO: handle error on ptarray construction */
+               /* XXX jorgearevalo: the error conditions aren't managed in ptarray_construct */
+               if (!rings[0]) {
+                       rterror("rt_raster_get_convex_hull: Unable to construct point array");
+                       return ES_ERROR;
+               }
+               pts = rings[0];
+
+               /* Upper-left corner (first and last points) */
+               p4d.x = gt[0];
+               p4d.y = gt[3];
+               ptarray_set_point4d(pts, 0, &p4d);
+               ptarray_set_point4d(pts, 4, &p4d);
+
+               /* Upper-right corner (we go clockwise) */
+               rt_raster_cell_to_geopoint(
+                       raster,
+                       raster->width, 0,
+                       &p4d.x, &p4d.y,
+                       gt
+               );
+               ptarray_set_point4d(pts, 1, &p4d);
+
+               /* Lower-right corner */
+               rt_raster_cell_to_geopoint(
+                       raster,
+                       raster->width, raster->height,
+                       &p4d.x, &p4d.y,
+                       gt
+               );
+               ptarray_set_point4d(pts, 2, &p4d);
+
+               /* Lower-left corner */
+               rt_raster_cell_to_geopoint(
+                       raster,
+                       0, raster->height,
+                       &p4d.x, &p4d.y,
+                       gt
+               );
+               ptarray_set_point4d(pts, 3, &p4d);
+
+               poly = lwpoly_construct(srid, 0, 1, rings);
+               *hull = lwpoly_as_lwgeom(poly);
+       }
+
+       return ES_NONE;
 }
 
 /**
@@ -6465,7 +6528,7 @@ rt_raster_compute_skewed_raster(
        int x;
        int y;
 
-       LWPOLY *spoly = NULL;
+       LWGEOM *geom = NULL;
        GEOSGeometry *sgeom = NULL;
        GEOSGeometry *ngeom = NULL;
 
@@ -6723,17 +6786,16 @@ rt_raster_compute_skewed_raster(
        do {
                covers = 0;
 
-               /* construct spoly from raster */
-               spoly = rt_raster_get_convex_hull(raster);
-               if (spoly == NULL) {
+               /* construct sgeom from raster */
+               if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
                        rterror("rt_raster_compute_skewed_raster: Unable to build skewed extent's geometry for covers test");
                        GEOSGeom_destroy(ngeom);
                        rt_raster_destroy(raster);
                        return NULL;
                }
 
-               sgeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(spoly));
-               lwpoly_free(spoly);
+               sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom);
+               lwgeom_free(geom);
 
                covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
                GEOSGeom_destroy(sgeom);
@@ -6772,17 +6834,16 @@ rt_raster_compute_skewed_raster(
                        else
                                raster->height--;
                        
-                       /* construct spoly from raster */
-                       spoly = rt_raster_get_convex_hull(raster);
-                       if (spoly == NULL) {
+                       /* construct sgeom from raster */
+                       if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
                                rterror("rt_raster_compute_skewed_raster: Unable to build skewed extent's geometry for minimizing dimensions");
                                GEOSGeom_destroy(ngeom);
                                rt_raster_destroy(raster);
                                return NULL;
                        }
 
-                       sgeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(spoly));
-                       lwpoly_free(spoly);
+                       sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom);
+                       lwgeom_free(geom);
 
                        covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
                        GEOSGeom_destroy(sgeom);
@@ -11615,7 +11676,7 @@ rt_raster_intersects(
        int j;
        int within = 0;
 
-       LWPOLY *hull[2] = {NULL};
+       LWGEOM *hull[2] = {NULL};
        GEOSGeometry *ghull[2] = {NULL};
 
        uint16_t width1;
@@ -11686,22 +11747,21 @@ rt_raster_intersects(
 
                rtn = 1;
                for (i = 0; i < 2; i++) {
-                       hull[i] = rt_raster_get_convex_hull(i < 1 ? rast1 : rast2);
-                       if (NULL == hull[i]) {
+                       if ((rt_raster_get_convex_hull(i < 1 ? rast1 : rast2, &(hull[i])) != ES_NONE) || NULL == hull[i]) {
                                for (j = 0; j < i; j++) {
                                        GEOSGeom_destroy(ghull[j]);
-                                       lwpoly_free(hull[j]);
+                                       lwgeom_free(hull[j]);
                                }
                                rtn = 0;
                                break;
                        }
-                       ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(hull[i]));
+                       ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(hull[i]);
                        if (NULL == ghull[i]) {
                                for (j = 0; j < i; j++) {
                                        GEOSGeom_destroy(ghull[j]);
-                                       lwpoly_free(hull[j]);
+                                       lwgeom_free(hull[j]);
                                }
-                               lwpoly_free(hull[i]);
+                               lwgeom_free(hull[i]);
                                rtn = 0;
                                break;
                        }
@@ -11722,7 +11782,7 @@ rt_raster_intersects(
 
                for (i = 0; i < 2; i++) {
                        GEOSGeom_destroy(ghull[i]);
-                       lwpoly_free(hull[i]);
+                       lwgeom_free(hull[i]);
                }
 
                if (rtn != 2) {
@@ -12912,9 +12972,8 @@ rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
        *surface = NULL;
 
        /* raster is empty, surface = NULL */
-       if (rt_raster_is_empty(raster)) {
+       if (rt_raster_is_empty(raster))
                return ES_NONE;
-       }
 
        /* if nband < 0, return the convex hull as a multipolygon */
        if (nband < 0) {
@@ -12924,7 +12983,10 @@ rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
                        hence the deep clone of the output geometry for returning
                        is the only way to guarentee the memory isn't shared
                */
-               tmp = lwpoly_as_lwgeom(rt_raster_get_convex_hull(raster));
+               if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
+                       rterror("rt_raster_surface: Unable to get convex hull of raster");
+                       return ES_ERROR;
+               }
                mpoly = lwgeom_as_multi(tmp);
                clone = lwgeom_clone_deep(mpoly);
                lwgeom_free(tmp);
@@ -12954,7 +13016,10 @@ rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
                        hence the deep clone of the output geometry for returning
                        is the only way to guarentee the memory isn't shared
                */
-               tmp = lwpoly_as_lwgeom(rt_raster_get_convex_hull(raster));
+               if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
+                       rterror("rt_raster_surface: Unable to get convex hull of raster");
+                       return ES_ERROR;
+               }
                mpoly = lwgeom_as_multi(tmp);
                clone = lwgeom_clone_deep(mpoly);
                lwgeom_free(tmp);

diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h
index 7c9fc4c486509e2cf9897f7f60da03b3e731f10b..a5578a3a83a5df6ecf01e4553af6eb959efe82e1 100644 (file)
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -1311,17 +1311,18 @@ rt_errorstate rt_raster_geopoint_to_cell(
 );
 
 /**
- * Get raster's polygon convex hull.
+ * Get raster's convex hull.
  *
- * The convex hull is a 4 vertices (5 to be closed) single
- * ring polygon bearing the raster's rotation
- * and using projection coordinates
+ * The convex hull is typically a 4 vertices (5 to be closed)
+ * single ring polygon bearing the raster's rotation and using
+ * projection coordinates.
  *
  * @param raster : the raster to get info from
+ * @param **hull : pointer to convex hull
  *
- * @return the convex hull, or NULL on error.
+ * @return ES_NONE if success, ES_ERROR if error
  */
-LWPOLY* rt_raster_get_convex_hull(rt_raster raster);
+rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull);
 
 /**
  * Get raster's envelope.

diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index f0da14e22ac25ebe12e63fed2a26590e3d50d525..50e7b9a061983d54a5e9dec2054a58194df2aae2 100644 (file)
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -943,51 +943,49 @@ Datum RASTER_to_binary(PG_FUNCTION_ARGS)
 }
 
 /**
- * Return the convex hull of this raster as a 4-vertices POLYGON.
+ * Return the convex hull of this raster
  */
 PG_FUNCTION_INFO_V1(RASTER_convex_hull);
 Datum RASTER_convex_hull(PG_FUNCTION_ARGS)
 {
-    rt_pgraster *pgraster;
-    rt_raster raster;
-    LWPOLY* convexhull = NULL;
-    GSERIALIZED* gser = NULL;
+       rt_pgraster *pgraster;
+       rt_raster raster;
+       LWGEOM *hull = NULL;
+       GSERIALIZED* gser = NULL;
+       size_t gser_size;
+       int err = ES_NONE;
 
-    if (PG_ARGISNULL(0)) PG_RETURN_NULL();
-    pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
+       if (PG_ARGISNULL(0))
+               PG_RETURN_NULL();
+       pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
 
-    {
-        raster = rt_raster_deserialize(pgraster, TRUE);
-        if ( ! raster ) {
-            elog(ERROR, "RASTER_convex_hull: Could not deserialize raster");
-                                               PG_FREE_IF_COPY(pgraster, 0);
-            PG_RETURN_NULL();
-        }
+       raster = rt_raster_deserialize(pgraster, TRUE);
+       if (!raster) {
+               elog(ERROR, "RASTER_convex_hull: Could not deserialize raster");
+               PG_FREE_IF_COPY(pgraster, 0);
+               PG_RETURN_NULL();
+       }
 
-        convexhull = rt_raster_get_convex_hull(raster);
-        if ( ! convexhull ) {
-            elog(ERROR, "RASTER_convex_hull: Could not get raster's convex hull");
-            rt_raster_destroy(raster);
-                                               PG_FREE_IF_COPY(pgraster, 0);
-            PG_RETURN_NULL();
-        }
-    }
+       err = rt_raster_get_convex_hull(raster, &hull);
+       rt_raster_destroy(raster);
+       PG_FREE_IF_COPY(pgraster, 0);
 
-    {
-        size_t gser_size;
-        gser = gserialized_from_lwgeom(lwpoly_as_lwgeom(convexhull), 0, &gser_size);
-        SET_VARSIZE(gser, gser_size);
-    }
+       if (err != ES_NONE) {
+               elog(ERROR, "RASTER_convex_hull: Could not get raster's convex hull");
+               PG_RETURN_NULL();
+       }
+       else if (hull == NULL) {
+               elog(NOTICE, "Raster's convex hull is NULL");
+               PG_RETURN_NULL();
+       }
 
-    /* Free raster and lwgeom memory */
-    rt_raster_destroy(raster);
-    pfree(convexhull);
-               PG_FREE_IF_COPY(pgraster, 0);
+       gser = gserialized_from_lwgeom(hull, 0, &gser_size);
+       lwgeom_free(hull);
 
-    PG_RETURN_POINTER(gser);
+       SET_VARSIZE(gser, gser_size);
+       PG_RETURN_POINTER(gser);
 }
 
-
 PG_FUNCTION_INFO_V1(RASTER_dumpAsPolygons);
 Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS) {
        FuncCallContext *funcctx;
@@ -16784,8 +16782,7 @@ Datum RASTER_clip(PG_FUNCTION_ARGS)
                arg->extenttype = ET_FIRST;
 
        /* get intersection geometry of input raster and input geometry */
-       rastgeom = lwpoly_as_lwgeom(rt_raster_get_convex_hull(arg->raster));
-       if (rastgeom == NULL) {
+       if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
                elog(ERROR, "RASTER_clip: Unable to get convex hull of raster");
 
                rtpg_clip_arg_destroy(arg);

diff --git a/raster/test/cunit/cu_raster_geometry.c b/raster/test/cunit/cu_raster_geometry.c
index 7d9ff15f8c813c974e8142d2108973c36bfde851..1f50015c5750cba36d27dd074371a172afbdb55a 100644 (file)
--- a/raster/test/cunit/cu_raster_geometry.c
+++ b/raster/test/cunit/cu_raster_geometry.c
@@ -26,11 +26,40 @@
 
 static void test_raster_convex_hull() {
        rt_raster raster = NULL;
-       LWPOLY *convexhull = NULL;
+       LWGEOM *hull = NULL;
+       LWPOLY *poly = NULL;
        POINTARRAY *ring = NULL;
        POINT4D pt;
 
-       /* create raster */
+       /* NULL raster */
+       CU_ASSERT_EQUAL(rt_raster_get_convex_hull(NULL, &hull), ES_NONE);
+       CU_ASSERT(hull == NULL);
+
+       /* width = 0, height = 0 */
+       raster = rt_raster_new(0, 0);
+       CU_ASSERT(raster != NULL);
+
+       CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
+       CU_ASSERT_EQUAL(hull->type, POINTTYPE);
+       cu_free_raster(raster);
+
+       /* width = 0 */
+       raster = rt_raster_new(0, 256);
+       CU_ASSERT(raster != NULL);
+
+       CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
+       CU_ASSERT_EQUAL(hull->type, LINETYPE);
+       cu_free_raster(raster);
+
+       /* height = 0 */
+       raster = rt_raster_new(256, 0);
+       CU_ASSERT(raster != NULL);
+
+       CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
+       CU_ASSERT_EQUAL(hull->type, LINETYPE);
+       cu_free_raster(raster);
+
+       /* normal raster */
        raster = rt_raster_new(256, 256);
        CU_ASSERT(raster != NULL);
 
@@ -38,11 +67,12 @@ static void test_raster_convex_hull() {
        rt_raster_set_scale(raster, 1, 1);
        rt_raster_set_skews(raster, 4, 5);
 
-       convexhull = rt_raster_get_convex_hull(raster);
-       CU_ASSERT_EQUAL(convexhull->srid, rt_raster_get_srid(raster));
-       CU_ASSERT_EQUAL(convexhull->nrings, 1);
+       CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
+       poly = lwgeom_as_lwpoly(hull);
+       CU_ASSERT_EQUAL(poly->srid, rt_raster_get_srid(raster));
+       CU_ASSERT_EQUAL(poly->nrings, 1);
 
-       ring = convexhull->rings[0];
+       ring = poly->rings[0];
        CU_ASSERT(ring != NULL);
        CU_ASSERT_EQUAL(ring->npoints, 5);
 
@@ -66,7 +96,7 @@ static void test_raster_convex_hull() {
        CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
        CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
 
-       lwpoly_free(convexhull);
+       lwgeom_free(hull);
        cu_free_raster(raster);
 }