From cebb5c8923bd4a619da679b5e2e5bc4f451a39cb Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Fri, 7 Dec 2012 01:08:21 +0000 Subject: [PATCH] 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 --- NEWS | 1 + raster/rt_core/rt_api.c | 227 ++++++++++++++++--------- raster/rt_core/rt_api.h | 13 +- raster/rt_pg/rt_pg.c | 67 ++++---- raster/test/cunit/cu_raster_geometry.c | 44 ++++- 5 files changed, 223 insertions(+), 129 deletions(-) diff --git a/NEWS b/NEWS index 6652b19c6..1482a756f 100644 --- 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 6ef2f3774..06a3a2a3f 100644 --- 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 7c9fc4c48..a5578a3a8 100644 --- 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 f0da14e22..50e7b9a06 100644 --- 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 7d9ff15f8..1f50015c5 100644 --- 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); } -- 2.50.1