From: Bborie Park Date: Tue, 26 Feb 2013 23:23:38 +0000 (+0000) Subject: Added ST_MinConvexHull(raster). Ticket #2210 X-Git-Tag: 2.1.0beta2~191 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=eb75a33623594b8b4c5b584d844c6f02464ed637;p=postgis Added ST_MinConvexHull(raster). Ticket #2210 git-svn-id: http://svn.osgeo.org/postgis/trunk@11118 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index b65a63ea0..44da35448 100644 --- a/NEWS +++ b/NEWS @@ -67,6 +67,7 @@ PostGIS 2.1.0 - #1293, ST_Resize(raster) to resize rasters based upon width/height - #945, improved join selectivity, N-D selectivity calculations, user accessible selectivity and stats reader functions for testing + - #2210, ST_MinConvexHull(raster) * Enhancements * - #823, tiger geocoder: Make loader_generate_script download portion diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index c042b9e48..46105d951 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -7731,30 +7731,17 @@ WHERE rid = 4; Examples - Refer to PostGIS Raster Specification - for a diagram of this. - - - - - - --- Note envelope and convexhull are more or less the same + Refer to PostGIS Raster Specification for a diagram of this. + +-- Note envelope and convexhull are more or less the same SELECT ST_AsText(ST_ConvexHull(rast)) As convhull, ST_AsText(ST_Envelope(rast)) As env FROM dummy_rast WHERE rid=1; convhull | env - ---------------------------------------------------------+----------------------- +--------------------------------------------------------+------------------------------------ POLYGON((0.5 0.5,20.5 0.5,20.5 60.5,0.5 60.5,0.5 0.5)) | POLYGON((0 0,20 0,20 60,0 60,0 0)) - - - - - - - + -- now we skew the raster -- note how the convex hull and envelope are now different @@ -7764,25 +7751,109 @@ FROM (SELECT ST_SetRotation(rast, 0.1, 0.1) As rast FROM dummy_rast WHERE rid=1) As foo; convhull | env - --------------------------------------------------------+------------------------------------ POLYGON((0.5 0.5,20.5 1.5,22.5 61.5,2.5 60.5,0.5 0.5)) | POLYGON((0 0,22 0,22 61,0 61,0 0)) - - - - - - See Also - , , + + , + , + , + + + + + + ST_MinConvexHull + + Return the convex hull geometry of the raster excluding NODATA pixels. + + + + + + geometry ST_MinConvexHull + raster rast + integer nband=NULL + + + + + + Description + + + Return the convex hull geometry of the raster excluding NODATA pixels. If nband is NULL, all bands of the raster are considered. + + + Availability: 2.1.0 + + + + Examples + +WITH foo AS ( + SELECT + ST_SetValues( + ST_SetValues( + ST_AddBand(ST_AddBand(ST_MakeEmptyRaster(9, 9, 0, 0, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0), 2, '8BUI', 1, 0), + 1, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 1], + [0, 0, 0, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ), + 2, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ) AS rast +) +SELECT + ST_AsText(ST_ConvexHull(rast)) AS hull, + ST_AsText(ST_MinConvexHull(rast)) AS mhull, + ST_AsText(ST_MinConvexHull(rast, 1)) AS mhull_1, + ST_AsText(ST_MinConvexHull(rast, 2)) AS mhull_2 +FROM foo + + hull | mhull | mhull_1 | mhull_2 +----------------------------------+-------------------------------------+-------------------------------------+------------------------------------- + POLYGON((0 0,9 0,9 -9,0 -9,0 0)) | POLYGON((0 -3,9 -3,9 -9,0 -9,0 -3)) | POLYGON((3 -3,9 -3,9 -6,3 -6,3 -3)) | POLYGON((0 -3,6 -3,6 -9,0 -9,0 -3)) + + + + See Also + + , + , + , + + + + + ST_DumpAsPolygons diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 6ad24bbe3..8f124ff07 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -4,7 +4,7 @@ * WKTRaster - Raster Types for PostGIS * http://trac.osgeo.org/postgis/wiki/WKTRaster * - * Copyright (C) 2011-2012 Regents of the University of California + * Copyright (C) 2011-2013 Regents of the University of California * * Copyright (C) 2010-2011 Jorge Arevalo * Copyright (C) 2010-2011 David Zwarg @@ -14269,3 +14269,306 @@ rt_raster_iterator( *rtnraster = rtnrast; return ES_NONE; } + +/****************************************************************************** +* rt_raster_perimeter() +******************************************************************************/ +static rt_errorstate +_rti_raster_get_band_perimeter(rt_band band, uint16_t *trim) { + uint16_t width = 0; + uint16_t height = 0; + int x = 0; + int y = 0; + int offset = 0; + int done[4] = {0}; + double value = 0; + int nodata = 0; + + assert(band != NULL); + assert(band->raster != NULL); + assert(trim != NULL); + + memset(trim, 0, sizeof(uint16_t) * 4); + + width = rt_band_get_width(band); + height = rt_band_get_height(band); + + /* top */ + for (y = 0; y < height; y++) { + for (offset = 0; offset < 3; offset++) { + /* every third pixel */ + for (x = offset; x < width; x += 3) { + if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) { + rterror("_rti_raster_get_band_perimeter: Unable to get band pixel"); + return ES_ERROR; + } + + RASTER_DEBUGF(4, "top (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata); + if (!nodata) { + trim[0] = y; + done[0] = 1; + break; + } + } + + if (done[0]) + break; + } + + if (done[0]) + break; + } + + /* right */ + for (x = width - 1; x >= 0; x--) { + for (offset = 0; offset < 3; offset++) { + /* every third pixel */ + for (y = offset; y < height; y += 3) { + if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) { + rterror("_rti_raster_get_band_perimeter: Unable to get band pixel"); + return ES_ERROR; + } + + RASTER_DEBUGF(4, "right (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata); + if (!nodata) { + trim[1] = width - (x + 1); + done[1] = 1; + break; + } + } + + if (done[1]) + break; + } + + if (done[1]) + break; + } + + /* bottom */ + for (y = height - 1; y >= 0; y--) { + for (offset = 0; offset < 3; offset++) { + /* every third pixel */ + for (x = offset; x < width; x += 3) { + if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) { + rterror("_rti_raster_get_band_perimeter: Unable to get band pixel"); + return ES_ERROR; + } + + RASTER_DEBUGF(4, "bottom (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata); + if (!nodata) { + trim[2] = height - (y + 1); + done[2] = 1; + break; + } + } + + if (done[2]) + break; + } + + if (done[2]) + break; + } + + /* left */ + for (x = 0; x < width; x++) { + for (offset = 0; offset < 3; offset++) { + /* every third pixel */ + for (y = offset; y < height; y += 3) { + if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) { + rterror("_rti_raster_get_band_perimeter: Unable to get band pixel"); + return ES_ERROR; + } + + RASTER_DEBUGF(4, "left (x, , value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata); + if (!nodata) { + trim[3] = x; + done[3] = 1; + break; + } + } + + if (done[3]) + break; + } + + if (done[3]) + break; + } + + RASTER_DEBUGF(4, "trim = (%d, %d, %d, %d)", + trim[0], trim[1], trim[2], trim[3]); + + return ES_NONE; +} + +/** + * Get raster perimeter + * + * The perimeter is 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 nband : the band for the perimeter. 0-based + * value less than zero means all bands + * @param **perimeter : pointer to perimeter + * + * @return ES_NONE if success, ES_ERROR if error + */ +rt_errorstate rt_raster_get_perimeter( + rt_raster raster, int nband, + LWGEOM **perimeter +) { + rt_band band = NULL; + int numband = 0; + uint16_t *_nband = NULL; + int i = 0; + int j = 0; + uint16_t _trim[4] = {0}; + uint16_t trim[4] = {0}; /* top, right, bottom, left */ + int isset[4] = {0}; + double gt[6] = {0.0}; + int srid = SRID_UNKNOWN; + + POINTARRAY *pts = NULL; + POINT4D p4d; + POINTARRAY **rings = NULL; + LWPOLY* poly = NULL; + + assert(perimeter != NULL); + + *perimeter = NULL; + + /* empty raster, no perimeter */ + if (rt_raster_is_empty(raster)) + return ES_NONE; + + /* raster metadata */ + srid = rt_raster_get_srid(raster); + rt_raster_get_geotransform_matrix(raster, gt); + numband = rt_raster_get_num_bands(raster); + + RASTER_DEBUGF(3, "rt_raster_get_perimeter: raster is %dx%d", raster->width, raster->height); + + /* nband < 0 means all bands */ + if (nband >= 0) { + if (nband >= numband) { + rterror("rt_raster_get_boundary: Band %d not found for raster", nband); + return ES_ERROR; + } + + numband = 1; + } + else + nband = -1; + + RASTER_DEBUGF(3, "rt_raster_get_perimeter: nband, numband = %d, %d", nband, numband); + + _nband = rtalloc(sizeof(uint16_t) * numband); + if (_nband == NULL) { + rterror("rt_raster_get_boundary: Unable to allocate memory for band indices"); + return ES_ERROR; + } + + if (nband < 0) { + for (i = 0; i < numband; i++) + _nband[i] = i; + } + else + _nband[0] = nband; + + for (i = 0; i < numband; i++) { + band = rt_raster_get_band(raster, _nband[i]); + if (band == NULL) { + rterror("rt_raster_get_boundary: Unable to get band at index %d", _nband[i]); + rtdealloc(_nband); + return ES_ERROR; + } + + /* band is nodata */ + if (rt_band_get_isnodata_flag(band) != 0) + continue; + + if (_rti_raster_get_band_perimeter(band, trim) != ES_NONE) { + rterror("rt_raster_get_boundary: Unable band perimeter"); + rtdealloc(_nband); + return ES_ERROR; + } + + for (j = 0; j < 4; j++) { + if (!isset[j] || trim[j] < _trim[j]) { + _trim[j] = trim[j]; + isset[j] = 1; + } + } + } + + /* no longer needed */ + rtdealloc(_nband); + + /* check isset, just need to check one element */ + if (!isset[0]) { + /* return NULL as bands are empty */ + return ES_NONE; + } + + RASTER_DEBUGF(4, "trim = (%d, %d, %d, %d)", + trim[0], trim[1], trim[2], trim[3]); + + /* only one ring */ + rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*)); + if (!rings) { + rterror("rt_raster_get_perimeter: Unable to allocate memory for polygon ring"); + return ES_ERROR; + } + rings[0] = ptarray_construct(0, 0, 5); + if (!rings[0]) { + rterror("rt_raster_get_perimeter: Unable to construct point array"); + return ES_ERROR; + } + pts = rings[0]; + + /* Upper-left corner (first and last points) */ + rt_raster_cell_to_geopoint( + raster, + _trim[3], _trim[0], + &p4d.x, &p4d.y, + gt + ); + 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 - _trim[1], _trim[0], + &p4d.x, &p4d.y, + gt + ); + ptarray_set_point4d(pts, 1, &p4d); + + /* Lower-right corner */ + rt_raster_cell_to_geopoint( + raster, + raster->width - _trim[1], raster->height - _trim[2], + &p4d.x, &p4d.y, + gt + ); + ptarray_set_point4d(pts, 2, &p4d); + + /* Lower-left corner */ + rt_raster_cell_to_geopoint( + raster, + _trim[3], raster->height - _trim[2], + &p4d.x, &p4d.y, + gt + ); + ptarray_set_point4d(pts, 3, &p4d); + + poly = lwpoly_construct(srid, 0, 1, rings); + *perimeter = lwpoly_as_lwgeom(poly); + + return ES_NONE; +} diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 647b9bfc8..ed8a877dd 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -4,7 +4,7 @@ * WKTRaster - Raster Types for PostGIS * http://www.postgis.org/support/wiki/index.php?WKTRasterHomePage * - * Copyright (C) 2011-2012 Regents of the University of California + * Copyright (C) 2011-2013 Regents of the University of California * * Copyright (C) 2010-2011 Jorge Arevalo * Copyright (C) 2010-2011 David Zwarg @@ -1313,7 +1313,7 @@ rt_errorstate rt_raster_geopoint_to_cell( /** * Get raster's convex hull. * - * The convex hull is typically a 4 vertices (5 to be closed) + * The convex hull is a 4 vertices (5 to be closed) * single ring polygon bearing the raster's rotation and using * projection coordinates. * @@ -1339,6 +1339,24 @@ rt_errorstate rt_raster_get_envelope( rt_envelope *env ); +/** + * Get raster perimeter + * + * The perimeter is 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 nband : the band for the perimeter. 0-based + * @param **perimeter : pointer to perimeter + * + * @return ES_NONE if success, ES_ERROR if error + */ +rt_errorstate rt_raster_get_perimeter( + rt_raster raster, int nband, + LWGEOM **perimeter +); + /* * Compute skewed extent that covers unskewed extent. * diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 09fdfd8b5..7445578f3 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -383,6 +383,9 @@ Datum RASTER_union_finalfn(PG_FUNCTION_ARGS); /* raster clip */ Datum RASTER_clip(PG_FUNCTION_ARGS); +/* raster perimeter */ +Datum RASTER_perimeter(PG_FUNCTION_ARGS); + /* string replacement function taken from * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 */ @@ -967,23 +970,55 @@ Datum RASTER_convex_hull(PG_FUNCTION_ARGS) { rt_pgraster *pgraster; rt_raster raster; - LWGEOM *hull = NULL; + LWGEOM *geom = NULL; GSERIALIZED* gser = NULL; size_t gser_size; int err = ES_NONE; + bool minhull = FALSE; + 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); + /* # of args */ + if (PG_NARGS() > 1) + minhull = TRUE; + + if (!minhull) { + pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t)); + raster = rt_raster_deserialize(pgraster, TRUE); + } + else { + pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + raster = rt_raster_deserialize(pgraster, FALSE); + } + if (!raster) { elog(ERROR, "RASTER_convex_hull: Could not deserialize raster"); PG_FREE_IF_COPY(pgraster, 0); PG_RETURN_NULL(); } - err = rt_raster_get_convex_hull(raster, &hull); + if (!minhull) + err = rt_raster_get_convex_hull(raster, &geom); + else { + int nband = -1; + + /* get arg 1 */ + if (!PG_ARGISNULL(1)) { + nband = PG_GETARG_INT32(1); + if (!rt_raster_has_band(raster, nband - 1)) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + nband = nband - 1; + } + + err = rt_raster_get_perimeter(raster, nband, &geom); + } + rt_raster_destroy(raster); PG_FREE_IF_COPY(pgraster, 0); @@ -991,13 +1026,13 @@ Datum RASTER_convex_hull(PG_FUNCTION_ARGS) elog(ERROR, "RASTER_convex_hull: Could not get raster's convex hull"); PG_RETURN_NULL(); } - else if (hull == NULL) { + else if (geom == NULL) { elog(NOTICE, "Raster's convex hull is NULL"); PG_RETURN_NULL(); } - gser = gserialized_from_lwgeom(hull, 0, &gser_size); - lwgeom_free(hull); + gser = gserialized_from_lwgeom(geom, 0, &gser_size); + lwgeom_free(geom); SET_VARSIZE(gser, gser_size); PG_RETURN_POINTER(gser); @@ -17439,6 +17474,16 @@ Datum RASTER_clip(PG_FUNCTION_ARGS) PG_RETURN_POINTER(pgrtn); } +/* ---------------------------------------------------------------- */ +/* Find perimeter of raster */ +/* ---------------------------------------------------------------- */ + +PG_FUNCTION_INFO_V1(RASTER_perimeter); +Datum RASTER_perimeter(PG_FUNCTION_ARGS) +{ + PG_RETURN_NULL(); +} + /* ---------------------------------------------------------------- */ /* Memory allocation / error reporting hooks */ /* TODO: reuse the ones in libpgcommon ? */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index e54b0f173..19d5f4417 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -103,6 +103,14 @@ CREATE OR REPLACE FUNCTION st_convexhull(raster) AS 'MODULE_PATHNAME','RASTER_convex_hull' LANGUAGE 'c' IMMUTABLE STRICT; +CREATE OR REPLACE FUNCTION st_minconvexhull( + rast raster, + nband integer DEFAULT NULL +) + RETURNS geometry + AS 'MODULE_PATHNAME','RASTER_convex_hull' + LANGUAGE 'c' IMMUTABLE; + CREATE OR REPLACE FUNCTION box3d(raster) RETURNS box3d AS 'select box3d(st_convexhull($1))' diff --git a/raster/test/cunit/cu_raster_geometry.c b/raster/test/cunit/cu_raster_geometry.c index 9569e28b8..c698e1416 100644 --- a/raster/test/cunit/cu_raster_geometry.c +++ b/raster/test/cunit/cu_raster_geometry.c @@ -2,7 +2,7 @@ * PostGIS Raster - Raster Types for PostGIS * http://www.postgis.org/support/wiki/index.php?WKTRasterHomePage * - * Copyright (C) 2012 Regents of the University of California + * Copyright (C) 2012-2013 Regents of the University of California * * * This program is free software; you can redistribute it and/or modify @@ -233,6 +233,160 @@ static void test_raster_surface() { cu_free_raster(rast); } +static void test_raster_perimeter() { + rt_raster rast; + rt_band band; + const int maxX = 5; + const int maxY = 5; + int x, y; + char *wkt = NULL; + LWGEOM *geom = NULL; + int err; + + rast = rt_raster_new(maxX, maxY); + CU_ASSERT(rast != NULL); + + rt_raster_set_offsets(rast, 0, 0); + rt_raster_set_scale(rast, 1, -1); + + band = cu_add_band(rast, PT_32BUI, 1, 0); + CU_ASSERT(band != NULL); + + for (y = 0; y < maxY; y++) { + for (x = 0; x < maxX; x++) { + rt_band_set_pixel(band, x, y, 1, NULL); + } + } + + err = rt_raster_get_perimeter(rast, -1, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,5 0,5 -5,0 -5,0 0))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* row 0 is NODATA */ + rt_band_set_pixel(band, 0, 0, 0, NULL); + rt_band_set_pixel(band, 1, 0, 0, NULL); + rt_band_set_pixel(band, 2, 0, 0, NULL); + rt_band_set_pixel(band, 3, 0, 0, NULL); + rt_band_set_pixel(band, 4, 0, 0, NULL); + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,5 -1,5 -5,0 -5,0 -1))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* column 4 is NODATA */ + /* pixel 4, 0 already set to NODATA */ + rt_band_set_pixel(band, 4, 1, 0, NULL); + rt_band_set_pixel(band, 4, 2, 0, NULL); + rt_band_set_pixel(band, 4, 3, 0, NULL); + rt_band_set_pixel(band, 4, 4, 0, NULL); + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,4 -1,4 -5,0 -5,0 -1))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* row 4 is NODATA */ + rt_band_set_pixel(band, 0, 4, 0, NULL); + rt_band_set_pixel(band, 1, 4, 0, NULL); + rt_band_set_pixel(band, 2, 4, 0, NULL); + rt_band_set_pixel(band, 3, 4, 0, NULL); + /* pixel 4, 4 already set to NODATA */ + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,4 -1,4 -4,0 -4,0 -1))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* column 0 is NODATA */ + /* pixel 0, 0 already set to NODATA*/ + rt_band_set_pixel(band, 0, 1, 0, NULL); + rt_band_set_pixel(band, 0, 2, 0, NULL); + rt_band_set_pixel(band, 0, 3, 0, NULL); + /* pixel 0, 4 already set to NODATA*/ + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((1 -1,4 -1,4 -4,1 -4,1 -1))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* columns 1 and 3 are NODATA */ + /* pixel 1, 0 already set to NODATA */ + rt_band_set_pixel(band, 1, 1, 0, NULL); + rt_band_set_pixel(band, 1, 2, 0, NULL); + rt_band_set_pixel(band, 1, 3, 0, NULL); + /* pixel 1, 4 already set to NODATA */ + /* pixel 3, 0 already set to NODATA */ + rt_band_set_pixel(band, 3, 1, 0, NULL); + rt_band_set_pixel(band, 3, 2, 0, NULL); + rt_band_set_pixel(band, 3, 3, 0, NULL); + /* pixel 3, 4 already set to NODATA */ + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((2 -1,3 -1,3 -4,2 -4,2 -1))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* more pixels are NODATA */ + rt_band_set_pixel(band, 2, 1, 0, NULL); + rt_band_set_pixel(band, 2, 3, 0, NULL); + + err = rt_raster_get_perimeter(rast, 0, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((2 -2,3 -2,3 -3,2 -3,2 -2))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + /* second band */ + band = cu_add_band(rast, PT_32BUI, 1, 0); + CU_ASSERT(band != NULL); + + for (y = 0; y < maxY; y++) { + for (x = 0; x < maxX; x++) { + rt_band_set_pixel(band, x, y, 1, NULL); + } + } + + err = rt_raster_get_perimeter(rast, -1, &geom); + CU_ASSERT_EQUAL(err, ES_NONE); + CU_ASSERT(geom != NULL); + wkt = lwgeom_to_text(geom); + CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,5 0,5 -5,0 -5,0 0))"); + rtdealloc(wkt); + lwgeom_free(geom); + geom = NULL; + + cu_free_raster(rast); +} + static void test_raster_pixel_as_polygon() { rt_raster rast; rt_band band; @@ -286,6 +440,7 @@ CU_TestInfo raster_geometry_tests[] = { */ PG_TEST(test_raster_convex_hull), PG_TEST(test_raster_surface), + PG_TEST(test_raster_perimeter), PG_TEST(test_raster_pixel_as_polygon), CU_TEST_INFO_NULL }; diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index d64d9c600..b0a2dac0d 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -3,7 +3,7 @@ # # Copyright (c) 2009 Sandro Santilli , Pierre Racine # Copyright (c) 2011 Jorge Arevalo -# Copyright (c) 2011-2012 Regents of the University of California +# Copyright (c) 2011-2013 Regents of the University of California # # # This is free software; you can redistribute and/or modify it under @@ -58,7 +58,8 @@ TEST_PROPS = \ rt_hasnoband \ rt_metadata \ rt_rastertoworldcoord \ - rt_worldtorastercoord + rt_worldtorastercoord \ + rt_convexhull TEST_BANDPROPS = \ rt_band_properties \ diff --git a/raster/test/regress/rt_convexhull.sql b/raster/test/regress/rt_convexhull.sql new file mode 100644 index 000000000..18496cd75 --- /dev/null +++ b/raster/test/regress/rt_convexhull.sql @@ -0,0 +1,136 @@ +DROP TABLE IF EXISTS test_raster_convexhull; + +CREATE TABLE test_raster_convexhull AS + SELECT ST_AddBand(ST_AddBand(ST_MakeEmptyRaster(9, 9, 0, 0, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0), 2, '8BUI', 1, 0) AS rast; + +SELECT + ST_AsText(ST_ConvexHull(rast)), + ST_AsText(ST_MinConvexHull(rast)), + ST_AsText(ST_MinConvexHull(rast, 1)), + ST_AsText(ST_MinConvexHull(rast, 2)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 1, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 1)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 1, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 1)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 1, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 1)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 1, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 1], + [0, 0, 0, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 1)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 2, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 2)), + ST_AsText(ST_MinConvexHull(rast)) +FROM test_raster_convexhull; + +UPDATE test_raster_convexhull SET + rast = ST_SetValues( + rast, 2, 1, 1, + ARRAY[ + [0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0] + ]::double precision[][] + ); + +SELECT + ST_AsText(ST_MinConvexHull(rast, 1)), + ST_AsText(ST_MinConvexHull(rast, 2)), + ST_AsText(ST_MinConvexHull(rast)) +FROM test_raster_convexhull; + +DROP TABLE IF EXISTS test_raster_convexhull; diff --git a/raster/test/regress/rt_convexhull_expected b/raster/test/regress/rt_convexhull_expected new file mode 100644 index 000000000..2142d6175 --- /dev/null +++ b/raster/test/regress/rt_convexhull_expected @@ -0,0 +1,9 @@ +NOTICE: table "test_raster_convexhull" does not exist, skipping +NOTICE: Raster's convex hull is NULL +POLYGON((0 0,9 0,9 -9,0 -9,0 0))|POLYGON((0 0,9 0,9 -9,0 -9,0 0))||POLYGON((0 0,9 0,9 -9,0 -9,0 0)) +POLYGON((4 -4,5 -4,5 -5,4 -5,4 -4)) +POLYGON((3 -3,5 -3,5 -6,3 -6,3 -3)) +POLYGON((3 -3,6 -3,6 -6,3 -6,3 -3)) +POLYGON((3 -3,9 -3,9 -6,3 -6,3 -3)) +POLYGON((0 -3,6 -3,6 -6,0 -6,0 -3))|POLYGON((0 -3,9 -3,9 -6,0 -6,0 -3)) +POLYGON((3 -3,9 -3,9 -6,3 -6,3 -3))|POLYGON((0 0,9 0,9 -9,0 -9,0 0))|POLYGON((0 0,9 0,9 -9,0 -9,0 0))