From: Bborie Park Date: Tue, 22 May 2012 17:05:14 +0000 (+0000) Subject: addition of ST_NearestValue and ST_Neighborhood X-Git-Tag: 2.1.0beta2~997 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=bc14c40e50884ac053639d7738854a2cb3435da0;p=postgis addition of ST_NearestValue and ST_Neighborhood git-svn-id: http://svn.osgeo.org/postgis/trunk@9786 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 9f9285b1d..0f9d3051a 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -2102,6 +2102,238 @@ rt_band_get_pixel( } } +/** + * Get nearest pixel(s) with value (not NODATA) to specified pixel + * + * @param band: the band to get nearest pixel(s) from + * @param x: the column of the pixel (0-based) + * @param y: the line of the pixel (0-based) + * @param distance: the number of pixels around the specified pixel + * @param exclude_nodata_value: if non-zero, ignore nodata values + * to check for pixels with value + * @param npixels: return set of rt_pixel object or NULL + * + * @return -1 on error, otherwise the number of rt_pixel objects + * in npixels + */ +int rt_band_get_nearest_pixel( + rt_band band, + int x, int y, + uint16_t distance, + int exclude_nodata_value, + rt_pixel *npixels +) { + rt_pixel npixel = NULL; + int extent[4] = {0}; + int d0 = 0; + uint32_t _distance; + uint32_t i = 0; + uint32_t j = 0; + uint32_t k = 0; + uint32_t max = 0; + int _x = 0; + int _y = 0; + int *_inc = NULL; + double pixval = 0; + double minval = 0; + int err; + uint32_t count = 0; + + assert(NULL != band); + assert(NULL != npixels); + + RASTER_DEBUG(3, "Starting"); + + minval = rt_pixtype_get_min_value(band->pixtype); + + if (!distance) + d0 = 1; + + RASTER_DEBUGF(4, "Selected pixel: %d x %d", x, y); + RASTER_DEBUGF(4, "Distance: %d", distance); + + /* shortcuts if outside band extent */ + if ( + exclude_nodata_value && ( + (x < 0 || x > band->width) || + (y < 0 || y > band->height) + ) + ) { + /* no distance specified, jump to pixel close to extent */ + if (!distance) { + if (x < 0) + x = -1; + else if (x > band->width) + x = band->width; + + if (y < 0) + y = -1; + else if (y > band->height) + y = band->height; + + RASTER_DEBUGF(4, "Moved selected pixel: %d x %d", x, y); + } + /* + distance specified + if distance won't capture extent of band, return 0 + */ + else if ( + ((x < 0 && abs(x) > distance) || (x - band->width >= distance)) || + ((x < 0 && abs(y) > distance) || (y - band->height >= distance)) + ) { + RASTER_DEBUG(4, "No nearest pixels possible for provided pixel and distance"); + return 0; + } + } + + /* determine the maximum distance to prevent an infinite loop */ + if (!distance) { + int a, b; + a = abs(x); + b = abs(x - band->width); + + if (a > b) + distance = a; + else + distance = b; + + a = abs(y); + b = abs(y - band->height); + if (a > distance) + distance = a; + if (b > distance) + distance = b; + + RASTER_DEBUGF(4, "Maximum distance: %d", distance); + } + + count = 0; + _distance = 0; + *npixels = NULL; + do { + _distance++; + extent[0] = x - _distance; /* min x */ + extent[1] = y - _distance; /* min y */ + extent[2] = x + _distance; /* max x */ + extent[3] = y + _distance; /* max y */ + + RASTER_DEBUGF(4, "Processing distance: %d", _distance); + RASTER_DEBUGF(4, "Extent: (%d, %d, %d, %d)", + extent[0], extent[1], extent[2], extent[3]); + + for (i = 0; i < 2; i++) { + + /* by row */ + if (i < 1) + max = extent[2] - extent[0] + 1; + /* by column */ + else + max = extent[3] - extent[1] + 1; + max = abs(max); + + for (j = 0; j < 2; j++) { + /* by row */ + if (i < 1) { + _x = extent[0]; + _inc = &_x; + + /* top row */ + if (j < 1) + _y = extent[1]; + /* bottom row */ + else + _y = extent[3]; + } + /* by column */ + else { + _y = extent[1] + 1; + _inc = &_y; + + /* left column */ + if (j < 1) { + _x = extent[0]; + max -= 2; + } + /* right column */ + else + _x = extent[2]; + } + + RASTER_DEBUGF(4, "max: %d", max); + for (k = 0; k < max; k++) { + /* outside band extent, set to NODATA */ + if ( + (_x < 0 || _x >= band->width) || + (_y < 0 || _y >= band->height) + ) { + /* no NODATA, set to minimum possible value */ + if (!band->hasnodata) + pixval = minval; + else + pixval = band->nodataval; + RASTER_DEBUGF(4, "NODATA pixel outside band extent: (x, y, val) = (%d, %d, %f)", _x, _y, pixval); + } + else { + err = rt_band_get_pixel( + band, + _x, _y, + &pixval + ); + if (err < 0) { + rterror("rt_band_get_nearest_pixel: Unable to get pixel value"); + if (count) rtdealloc(*npixels); + return -1; + } + RASTER_DEBUGF(4, "Pixel: (x, y, val) = (%d, %d, %f)", _x, _y, pixval); + } + + /* use pixval? */ + if ( + !exclude_nodata_value || ( + exclude_nodata_value && + (band->hasnodata != FALSE) && ( + FLT_NEQ(pixval, band->nodataval) && + (rt_band_clamped_value_is_nodata(band, pixval) != 1) + ) + ) + ) { + /* add pixel to result set */ + RASTER_DEBUGF(4, "Adding pixel to set of nearest pixels: (x, y, val) = (%d, %d, %f)", _x, _y, pixval); + count++; + + if (*npixels == NULL) + *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count); + else + *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct rt_pixel_t) * count); + if (*npixels == NULL) { + rterror("rt_band_get_nearest_pixel: Unable to allocate memory for nearest pixel(s)"); + return -1; + } + + npixel = &((*npixels)[count - 1]); + npixel->x = _x; + npixel->y = _y; + npixel->value = pixval; + } + + (*_inc)++; + } + } + } + + /* distance threshhold met */ + if (_distance >= distance) + break; + else if (d0 && count) + break; + } + while (1); + + RASTER_DEBUGF(3, "Nearest pixels in return: %d", count); + + return count; +} + double rt_band_get_nodata(rt_band band) { diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index e14f03f4d..9e2182c39 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -133,6 +133,7 @@ */ typedef struct rt_raster_t* rt_raster; typedef struct rt_band_t* rt_band; +typedef struct rt_pixel_t* rt_pixel; typedef struct rt_geomval_t* rt_geomval; typedef struct rt_bandstats_t* rt_bandstats; typedef struct rt_histogram_t* rt_histogram; @@ -536,6 +537,28 @@ int rt_band_get_pixel( double *result ); +/** + * Get nearest pixel(s) with value (not NODATA) to specified pixel + * + * @param band: the band to get nearest pixel(s) from + * @param x: the column of the pixel (0-based) + * @param y: the line of the pixel (0-based) + * @param distance: the number of pixels around the specified pixel + * @param exclude_nodata_value: if non-zero, ignore nodata values + * to check for pixels with value + * @param npixels: return set of rt_pixel object or NULL + * + * @return -1 on error, otherwise the number of rt_pixel objects + * in npixels + */ +int rt_band_get_nearest_pixel( + rt_band band, + int x, int y, + uint16_t distance, + int exclude_nodata_value, + rt_pixel *npixels +); + /** * Returns the minimal possible value for the band according to the pixel type. * @param band: the band to get info from @@ -1637,6 +1660,12 @@ struct rt_band_t { }; +struct rt_pixel_t { + int x; /* column */ + int y; /* line */ + double value; +}; + /* polygon as LWPOLY with associated value */ struct rt_geomval_t { LWPOLY *geom; diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index d3b39bb3b..41d7cd93a 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -202,6 +202,12 @@ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS); /* Get pixel geographical shape */ Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS); +/* Get nearest value to a point */ +Datum RASTER_nearestValue(PG_FUNCTION_ARGS); + +/* Get the neighborhood around a pixel */ +Datum RASTER_neighborhood(PG_FUNCTION_ARGS); + /* Raster and band creation */ Datum RASTER_makeEmpty(PG_FUNCTION_ARGS); Datum RASTER_addband(PG_FUNCTION_ARGS); @@ -2359,6 +2365,368 @@ Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS) PG_RETURN_POINTER(gser); } +/** + * Return nearest value to a point + */ +PG_FUNCTION_INFO_V1(RASTER_nearestValue); +Datum RASTER_nearestValue(PG_FUNCTION_ARGS) +{ + rt_pgraster *pgraster = NULL; + rt_raster raster = NULL; + rt_band band = NULL; + int bandindex = 1; + int num_bands = 0; + GSERIALIZED *geom; + bool exclude_nodata_value = TRUE; + LWGEOM *lwgeom; + LWPOINT *point = NULL; + POINT2D p; + + double x; + double y; + int count; + rt_pixel npixels = NULL; + double value = 0; + int hasvalue = 0; + + if (PG_ARGISNULL(0)) + PG_RETURN_NULL(); + pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + raster = rt_raster_deserialize(pgraster, FALSE); + if (!raster) { + elog(ERROR, "RASTER_nearestValue: Could not deserialize raster"); + PG_RETURN_NULL(); + } + + /* band index is 1-based */ + if (!PG_ARGISNULL(1)) + bandindex = PG_GETARG_INT32(1); + num_bands = rt_raster_get_num_bands(raster); + if (bandindex < 1 || bandindex > num_bands) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* point */ + geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2)); + if (gserialized_get_type(geom) != POINTTYPE) { + elog(NOTICE, "Geometry provided must be a point"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* exclude_nodata_value flag */ + if (!PG_ARGISNULL(3)) + exclude_nodata_value = PG_GETARG_BOOL(3); + + /* get band */ + band = rt_raster_get_band(raster, bandindex - 1); + if (!band) { + elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* process geometry */ + lwgeom = lwgeom_from_gserialized(geom); + + if (lwgeom_is_empty(lwgeom)) { + elog(NOTICE, "Geometry provided cannot be empty"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + point = lwgeom_as_lwpoint(lwgeom); + getPoint2d_p(point->point, 0, &p); + + if (!rt_raster_geopoint_to_cell( + raster, + p.x, p.y, + &x, &y, + NULL + )) { + elog(ERROR, "RASTER_nearestValue: Unable to compute pixel coordinates from spatial coordinates"); + rt_raster_destroy(raster); + lwgeom_free(lwgeom); + PG_RETURN_NULL(); + } + + /* get value at point */ + if ( + (x >= 0 && x < rt_raster_get_width(raster)) && + (y >= 0 && y < rt_raster_get_height(raster)) + ) { + if (rt_band_get_pixel(band, x, y, &value) < 0) { + elog(ERROR, "RASTER_nearestValue: Unable to get pixel value for band at index %d", bandindex); + rt_raster_destroy(raster); + lwgeom_free(lwgeom); + PG_RETURN_NULL(); + } + + /* value at point, return value */ + if ( + !exclude_nodata_value || ( + exclude_nodata_value && + rt_band_get_hasnodata_flag(band) && ( + FLT_NEQ(value, rt_band_get_nodata(band)) && + (rt_band_clamped_value_is_nodata(band, value) != 1) + ) + ) + ) { + rt_raster_destroy(raster); + lwgeom_free(lwgeom); + + PG_RETURN_FLOAT8(value); + } + } + + /* get neighborhood */ + count = rt_band_get_nearest_pixel( + band, + x, y, + 0, + exclude_nodata_value, + &npixels + ); + rt_band_destroy(band); + /* error or no neighbors */ + if (count < 1) { + /* error */ + if (count < 0) + elog(NOTICE, "Unable to get the nearest value for band at index %d", bandindex); + /* no nearest pixel */ + else + elog(NOTICE, "No nearest value found for band at index %d", bandindex); + + lwgeom_free(lwgeom); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* more than one nearest value, see which one is closest */ + if (count > 1) { + int i = 0; + LWPOLY *poly = NULL; + double lastdist = -1; + double dist; + + for (i = 0; i < count; i++) { + /* convex-hull of pixel */ + poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y); + if (!poly) { + elog(ERROR, "RASTER_nearestValue: Unable to get polygon of neighboring pixel"); + lwgeom_free(lwgeom); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* distance between convex-hull and point */ + dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom); + if (lastdist < 0 || dist < lastdist) { + value = npixels[i].value; + hasvalue = 1; + } + lastdist = dist; + + lwpoly_free(poly); + } + } + else { + value = npixels[0].value; + hasvalue = 1; + } + + pfree(npixels); + lwgeom_free(lwgeom); + rt_raster_destroy(raster); + + if (hasvalue) + PG_RETURN_FLOAT8(value); + else + PG_RETURN_NULL(); +} + +/** + * Return the neighborhood around a pixel + */ +PG_FUNCTION_INFO_V1(RASTER_neighborhood); +Datum RASTER_neighborhood(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + TupleDesc tupdesc; + + int call_cntr; + int max_calls; + + rt_pixel npixel1 = NULL; + rt_pixel npixel2 = NULL; + + if (SRF_IS_FIRSTCALL()) { + MemoryContext oldcontext; + + rt_pgraster *pgraster = NULL; + rt_raster raster = NULL; + rt_band band = NULL; + int bandindex = 1; + int num_bands = 0; + int x = 0; + int y = 0; + int distance = 0; + bool exclude_nodata_value = TRUE; + + int count; + + /* create a function context for cross-call persistence */ + funcctx = SRF_FIRSTCALL_INIT(); + + /* switch to memory context appropriate for multiple function calls */ + oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + /* pgraster is null, return nothing */ + if (PG_ARGISNULL(0)) { + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + raster = rt_raster_deserialize(pgraster, FALSE); + if (!raster) { + elog(ERROR, "RASTER_neighborhood: Could not deserialize raster"); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + /* band index is 1-based */ + if (!PG_ARGISNULL(1)) + bandindex = PG_GETARG_INT32(1); + num_bands = rt_raster_get_num_bands(raster); + if (bandindex < 1 || bandindex > num_bands) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + /* pixel column, 1-based */ + x = PG_GETARG_INT32(2); + + /* pixel row, 1-based */ + y = PG_GETARG_INT32(3); + + /* distance */ + distance = PG_GETARG_INT32(4); + if (distance < 1) { + elog(NOTICE, "Invalid value for distance (must be greater than zero). Returning NULL"); + rt_raster_destroy(raster); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + /* exclude_nodata_value flag */ + if (!PG_ARGISNULL(5)) + exclude_nodata_value = PG_GETARG_BOOL(5); + + /* get band */ + band = rt_raster_get_band(raster, bandindex - 1); + if (!band) { + elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex); + rt_raster_destroy(raster); + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + /* get neighborhood */ + count = rt_band_get_nearest_pixel( + band, + x - 1, y - 1, + (uint16_t) distance, + exclude_nodata_value, + &npixel1 + ); + rt_band_destroy(band); + rt_raster_destroy(raster); + /* error or no neighbors */ + if (count < 1) { + /* error */ + if (count < 0) + elog(NOTICE, "Unable to get the pixel's neighborhood for band at index %d", bandindex); + /* no neighbors */ + else + elog(NOTICE, "Pixel has no neighbors for band at distance %d", distance); + + MemoryContextSwitchTo(oldcontext); + SRF_RETURN_DONE(funcctx); + } + + /* Store needed information */ + funcctx->user_fctx = npixel1; + + /* total number of tuples to be returned */ + funcctx->max_calls = count; + + /* Build a tuple descriptor for our result type */ + if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) { + ereport(ERROR, ( + errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg( + "function returning record called in context " + "that cannot accept type record" + ) + )); + } + + BlessTupleDesc(tupdesc); + funcctx->tuple_desc = tupdesc; + + MemoryContextSwitchTo(oldcontext); + } + + /* stuff done on every call of the function */ + funcctx = SRF_PERCALL_SETUP(); + + call_cntr = funcctx->call_cntr; + max_calls = funcctx->max_calls; + tupdesc = funcctx->tuple_desc; + npixel2 = funcctx->user_fctx; + + /* do when there is more left to send */ + if (call_cntr < max_calls) { + int values_length = 3; + Datum values[values_length]; + bool *nulls = NULL; + HeapTuple tuple; + Datum result; + + POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr); + + nulls = palloc(sizeof(bool) * values_length); + memset(nulls, FALSE, values_length); + + /* x,y are 0-based, make 1-based for end users */ + values[0] = Int64GetDatum(npixel2[call_cntr].x + 1); + values[1] = Int64GetDatum(npixel2[call_cntr].y + 1); + values[2] = Float8GetDatum(npixel2[call_cntr].value); + + /* build a tuple */ + tuple = heap_form_tuple(tupdesc, values, nulls); + + /* make the tuple into a datum */ + result = HeapTupleGetDatum(tuple); + + /* clean up */ + pfree(nulls); + + SRF_RETURN_NEXT(funcctx, result); + } + /* do when there is no more left */ + else { + pfree(npixel2); + SRF_RETURN_DONE(funcctx); + } + +} + /** * Add a band to the given raster at the given position. */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index d1e44423a..f6e6846ed 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -3492,6 +3492,118 @@ CREATE OR REPLACE FUNCTION st_clip(rast raster, geom geometry, crop boolean) $$ SELECT ST_Clip($1, NULL, $2, null::float8[], $3) $$ LANGUAGE 'sql' STABLE; +----------------------------------------------------------------------- +-- ST_NearestValue +----------------------------------------------------------------------- + +CREATE OR REPLACE FUNCTION st_nearestvalue( + rast raster, band integer, + pt geometry, + exclude_nodata_value boolean DEFAULT TRUE +) + RETURNS double precision + AS 'MODULE_PATHNAME', 'RASTER_nearestValue' + LANGUAGE 'C' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_nearestvalue( + rast raster, + pt geometry, + exclude_nodata_value boolean DEFAULT TRUE +) + RETURNS double precision + AS $$ SELECT st_nearestvalue($1, 1, $2, $3) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_nearestvalue( + rast raster, band integer, + x integer, y integer, + exclude_nodata_value boolean DEFAULT TRUE +) + RETURNS double precision + AS $$ SELECT st_nearestvalue($1, $2, st_makepoint(st_raster2worldcoordx($1, $3, $4), st_raster2worldcoordy($1, $3, $4)), $5) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_nearestvalue( + rast raster, + x integer, y integer, + exclude_nodata_value boolean DEFAULT TRUE +) + RETURNS double precision + AS $$ SELECT st_nearestvalue($1, 1, st_makepoint(st_raster2worldcoordx($1, $2, $3), st_raster2worldcoordy($1, $2, $3)), $4) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +----------------------------------------------------------------------- +-- ST_Neighborhood +----------------------------------------------------------------------- + +CREATE OR REPLACE FUNCTION st_neighborhood( + rast raster, band integer, + ix integer, iy integer, + distance integer, + exclude_nodata_value boolean DEFAULT TRUE, + OUT x integer, OUT y integer, + OUT val double precision +) + RETURNS SETOF record + AS 'MODULE_PATHNAME', 'RASTER_neighborhood' + LANGUAGE 'C' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_neighborhood( + rast raster, + ix integer, iy integer, + distance integer, + exclude_nodata_value boolean DEFAULT TRUE, + OUT x integer, OUT y integer, + OUT val double precision +) + RETURNS SETOF record + AS $$ SELECT x, y, val FROM st_neighborhood($1, 1, $2, $3, $4, $5) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_neighborhood( + rast raster, band integer, + pt geometry, + distance integer, + exclude_nodata_value boolean DEFAULT TRUE, + OUT x integer, OUT y integer, + OUT val double precision +) + RETURNS SETOF record + AS $$ + DECLARE + wx int; + wy int; + BEGIN + IF (st_geometrytype($3) != 'ST_Point') THEN + RAISE EXCEPTION 'Attempting to get the neighbor of a pixel with a non-point geometry'; + END IF; + wx := st_x($3); + wy := st_y($3); + + RETURN QUERY + SELECT x, y, val + FROM st_neighborhood( + $1, $2, + st_world2rastercoordx(rast, wx, wy), + st_world2rastercoordy(rast, wx, wy), + $4, + $5 + ); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_neighborhood( + rast raster, + pt geometry, + distance integer, + exclude_nodata_value boolean DEFAULT TRUE, + OUT x integer, OUT y integer, + OUT val double precision +) + RETURNS SETOF record + AS $$ SELECT x, y, val FROM st_neighborhood($1, 1, $2, $3, $4) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + ------------------------------------------------------------------------------ -- raster constraint functions ------------------------------------------------------------------------------- diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index b2a29a7e4..1baebb177 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1225,9 +1225,8 @@ static void testBandStats() { raster = rt_raster_new(xmax, ymax); assert(raster); - band = addBand(raster, PT_32BUI, 0, 0); + band = addBand(raster, PT_32BUI, 1, 0); CHECK(band); - rt_band_set_nodata(band, 0); for (x = 0; x < xmax; x++) { for (y = 0; y < ymax; y++) { @@ -2590,6 +2589,193 @@ static void testCellGeoPoint() { deepRelease(raster); } +static void testNearestPixel() { + rt_raster rast; + rt_band band; + uint32_t x, y; + int rtn; + const int maxX = 10; + const int maxY = 10; + rt_pixel npixels = NULL; + + rast = rt_raster_new(maxX, maxY); + assert(rast); + + band = addBand(rast, PT_32BUI, 1, 0); + CHECK(band); + + for (x = 0; x < maxX; x++) { + for (y = 0; y < maxY; y++) { + rtn = rt_band_set_pixel(band, x, y, 1); + CHECK((rtn != -1)); + } + } + + rt_band_set_pixel(band, 0, 0, 0); + rt_band_set_pixel(band, 3, 0, 0); + rt_band_set_pixel(band, 6, 0, 0); + rt_band_set_pixel(band, 9, 0, 0); + rt_band_set_pixel(band, 1, 2, 0); + rt_band_set_pixel(band, 4, 2, 0); + rt_band_set_pixel(band, 7, 2, 0); + rt_band_set_pixel(band, 2, 4, 0); + rt_band_set_pixel(band, 5, 4, 0); + rt_band_set_pixel(band, 8, 4, 0); + rt_band_set_pixel(band, 0, 6, 0); + rt_band_set_pixel(band, 3, 6, 0); + rt_band_set_pixel(band, 6, 6, 0); + rt_band_set_pixel(band, 9, 6, 0); + rt_band_set_pixel(band, 1, 8, 0); + rt_band_set_pixel(band, 4, 8, 0); + rt_band_set_pixel(band, 7, 8, 0); + + /* 0,0 */ + rtn = rt_band_get_nearest_pixel( + band, + 0, 0, + 0, + 1, + &npixels + ); + CHECK((rtn == 3)); + if (rtn) + rtdealloc(npixels); + + /* 1,1 */ + rtn = rt_band_get_nearest_pixel( + band, + 1, 1, + 0, + 1, + &npixels + ); + CHECK((rtn == 6)); + if (rtn) + rtdealloc(npixels); + + /* 4,4 */ + rtn = rt_band_get_nearest_pixel( + band, + 4, 4, + 0, + 1, + &npixels + ); + CHECK((rtn == 7)); + if (rtn) + rtdealloc(npixels); + + /* 4,4 distance 2 */ + rtn = rt_band_get_nearest_pixel( + band, + 4, 4, + 2, + 1, + &npixels + ); + CHECK((rtn == 19)); + if (rtn) + rtdealloc(npixels); + + /* 10,10 */ + rtn = rt_band_get_nearest_pixel( + band, + 10, 10, + 0, + 1, + &npixels + ); + CHECK((rtn == 1)); + if (rtn) + rtdealloc(npixels); + + /* 11,11 distance 1 */ + rtn = rt_band_get_nearest_pixel( + band, + 11, 11, + 1, + 1, + &npixels + ); + CHECK((rtn == 0)); + if (rtn) + rtdealloc(npixels); + + /* -1,-1 */ + rtn = rt_band_get_nearest_pixel( + band, + -1, -1, + 0, + 1, + &npixels + ); + CHECK((rtn == 3)); + if (rtn) + rtdealloc(npixels); + + /* -1,-1 distance 1 */ + rtn = rt_band_get_nearest_pixel( + band, + -1, -1, + 1, + 1, + &npixels + ); + CHECK((rtn == 0)); + if (rtn) + rtdealloc(npixels); + + /* -1,1 distance 1 */ + rtn = rt_band_get_nearest_pixel( + band, + -1, 1, + 1, + 1, + &npixels + ); + CHECK((rtn == 2)); + if (rtn) + rtdealloc(npixels); + + /* -2,2 distance 1 */ + rtn = rt_band_get_nearest_pixel( + band, + -2, 2, + 1, + 1, + &npixels + ); + CHECK((rtn == 0)); + if (rtn) + rtdealloc(npixels); + + /* -10,2 distance 3 */ + rtn = rt_band_get_nearest_pixel( + band, + -10, 2, + 3, + 1, + &npixels + ); + CHECK((rtn == 0)); + if (rtn) + rtdealloc(npixels); + + /* -10,2 distance 3 include NODATA */ + rtn = rt_band_get_nearest_pixel( + band, + -10, 2, + 3, + 0, + &npixels + ); + CHECK((rtn == 48)); + if (rtn) + rtdealloc(npixels); + + deepRelease(rast); +} + int main() { @@ -2852,6 +3038,10 @@ main() testCellGeoPoint(); printf("OK\n"); + printf("Testing rt_band_get_nearest_pixel... "); + testNearestPixel(); + printf("OK\n"); + deepRelease(raster); return EXIT_SUCCESS; diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index e5df03980..230eb4448 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -74,7 +74,9 @@ TEST_BANDPROPS = \ rt_valuepercent \ rt_bandmetadata \ rt_pixelvalue \ - drop_rt_band_properties_test + drop_rt_band_properties_test \ + rt_neighborhood \ + rt_nearestvalue TEST_UTILITY = \ rt_utility \ diff --git a/raster/test/regress/rt_nearestvalue.sql b/raster/test/regress/rt_nearestvalue.sql new file mode 100644 index 000000000..43236ac0f --- /dev/null +++ b/raster/test/regress/rt_nearestvalue.sql @@ -0,0 +1,118 @@ +DROP TABLE IF EXISTS raster_nearestvalue; +CREATE TABLE raster_nearestvalue ( + rast raster +); +CREATE OR REPLACE FUNCTION make_test_raster() + RETURNS void + AS $$ + DECLARE + width int := 10; + height int := 10; + x int; + y int; + rast raster; + BEGIN + rast := ST_MakeEmptyRaster(width, height, 0, 0, 1, -1, 0, 0, 0); + rast := ST_AddBand(rast, 1, '8BUI', 1, 0); + + FOR y IN 1..height LOOP + FOR x IN 1..width LOOP + rast := ST_SetValue(rast, x, y, 2 * x + (1/3) * y); + END LOOP; + END LOOP; + + rast := ST_SetValue(rast, 1, 1, 0); + rast := ST_SetValue(rast, 4, 1, 0); + rast := ST_SetValue(rast, 7, 1, 0); + rast := ST_SetValue(rast, 10, 1, 0); + rast := ST_SetValue(rast, 2, 3, 0); + rast := ST_SetValue(rast, 5, 3, 0); + rast := ST_SetValue(rast, 8, 3, 0); + rast := ST_SetValue(rast, 3, 5, 0); + rast := ST_SetValue(rast, 6, 5, 0); + rast := ST_SetValue(rast, 9, 5, 0); + rast := ST_SetValue(rast, 1, 7, 0); + rast := ST_SetValue(rast, 4, 7, 0); + rast := ST_SetValue(rast, 7, 7, 0); + rast := ST_SetValue(rast, 10, 7, 0); + rast := ST_SetValue(rast, 2, 9, 0); + rast := ST_SetValue(rast, 5, 9, 0); + rast := ST_SetValue(rast, 8, 9, 0); + + INSERT INTO raster_nearestvalue VALUES (rast); + + RETURN; + END; + $$ LANGUAGE 'plpgsql'; +SELECT make_test_raster(); +DROP FUNCTION IF EXISTS make_test_raster(); + +SELECT + ST_NearestValue(rast, 1, 1, 1) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 2, 2) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 5, 5) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 5, 5) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 11, 11) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 12, 12) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 0, 0) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, 0, 2) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, -1, 3) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, -9, 3) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, -9, 3, FALSE) +FROM raster_nearestvalue; + +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(1, 1)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(2, 2)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(5, 5)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(5, 5)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(11, 11)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(12, 12)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(0, 0)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(0, 2)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(-1, 3)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(-9, 3)) +FROM raster_nearestvalue; +SELECT + ST_NearestValue(rast, 1, ST_MakePoint(-9, 3), FALSE) +FROM raster_nearestvalue; + +DROP TABLE IF EXISTS raster_nearestvalue; diff --git a/raster/test/regress/rt_nearestvalue_expected b/raster/test/regress/rt_nearestvalue_expected new file mode 100644 index 000000000..8ced9f4ee --- /dev/null +++ b/raster/test/regress/rt_nearestvalue_expected @@ -0,0 +1,23 @@ +NOTICE: table "raster_nearestvalue" does not exist, skipping +4 +4 +10 +10 +20 +20 +4 +2 +2 +2 +0 +4 +4 +10 +10 +18 +18 +4 +4 +4 +4 +0 diff --git a/raster/test/regress/rt_neighborhood.sql b/raster/test/regress/rt_neighborhood.sql new file mode 100644 index 000000000..1232dc9a3 --- /dev/null +++ b/raster/test/regress/rt_neighborhood.sql @@ -0,0 +1,74 @@ +DROP TABLE IF EXISTS raster_neighborhood; +CREATE TABLE raster_neighborhood ( + rast raster +); +CREATE OR REPLACE FUNCTION make_test_raster() + RETURNS void + AS $$ + DECLARE + rast raster; + BEGIN + rast := ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0); + rast := ST_AddBand(rast, 1, '8BUI', 1, 0); + + rast := ST_SetValue(rast, 1, 1, 0); + rast := ST_SetValue(rast, 4, 1, 0); + rast := ST_SetValue(rast, 7, 1, 0); + rast := ST_SetValue(rast, 10, 1, 0); + rast := ST_SetValue(rast, 2, 3, 0); + rast := ST_SetValue(rast, 5, 3, 0); + rast := ST_SetValue(rast, 8, 3, 0); + rast := ST_SetValue(rast, 3, 5, 0); + rast := ST_SetValue(rast, 6, 5, 0); + rast := ST_SetValue(rast, 9, 5, 0); + rast := ST_SetValue(rast, 1, 7, 0); + rast := ST_SetValue(rast, 4, 7, 0); + rast := ST_SetValue(rast, 7, 7, 0); + rast := ST_SetValue(rast, 10, 7, 0); + rast := ST_SetValue(rast, 2, 9, 0); + rast := ST_SetValue(rast, 5, 9, 0); + rast := ST_SetValue(rast, 8, 9, 0); + + INSERT INTO raster_neighborhood VALUES (rast); + + RETURN; + END; + $$ LANGUAGE 'plpgsql'; +SELECT make_test_raster(); +DROP FUNCTION IF EXISTS make_test_raster(); + +SELECT + ST_Neighborhood(rast, 1, 1, 1, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 2, 2, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 5, 5, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 5, 5, 2) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 11, 11, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 12, 12, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 0, 0, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, 0, 2, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, -1, 3, 1) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, -9, 3, 3) +FROM raster_neighborhood; +SELECT + ST_Neighborhood(rast, 1, -9, 3, 3, FALSE) +FROM raster_neighborhood; + +DROP TABLE IF EXISTS raster_neighborhood; diff --git a/raster/test/regress/rt_neighborhood_expected b/raster/test/regress/rt_neighborhood_expected new file mode 100644 index 000000000..0f0011d61 --- /dev/null +++ b/raster/test/regress/rt_neighborhood_expected @@ -0,0 +1,91 @@ +NOTICE: table "raster_neighborhood" does not exist, skipping +(1,2,1) +(2,2,1) +(2,1,1) +(2,1,1) +(3,1,1) +(1,3,1) +(3,3,1) +(1,2,1) +(3,2,1) +(4,4,1) +(5,4,1) +(6,4,1) +(4,6,1) +(5,6,1) +(6,6,1) +(4,5,1) +(4,4,1) +(5,4,1) +(6,4,1) +(4,6,1) +(5,6,1) +(6,6,1) +(4,5,1) +(3,3,1) +(4,3,1) +(6,3,1) +(7,3,1) +(3,7,1) +(5,7,1) +(6,7,1) +(3,4,1) +(3,6,1) +(7,4,1) +(7,5,1) +(7,6,1) +(10,10,1) +NOTICE: Pixel has no neighbors for band at distance 1 +NOTICE: Pixel has no neighbors for band at distance 1 +(1,3,1) +(1,2,1) +NOTICE: Pixel has no neighbors for band at distance 1 +NOTICE: Pixel has no neighbors for band at distance 3 +(-10,2,0) +(-9,2,0) +(-8,2,0) +(-10,4,0) +(-9,4,0) +(-8,4,0) +(-10,3,0) +(-8,3,0) +(-11,1,0) +(-10,1,0) +(-9,1,0) +(-8,1,0) +(-7,1,0) +(-11,5,0) +(-10,5,0) +(-9,5,0) +(-8,5,0) +(-7,5,0) +(-11,2,0) +(-11,3,0) +(-11,4,0) +(-7,2,0) +(-7,3,0) +(-7,4,0) +(-12,0,0) +(-11,0,0) +(-10,0,0) +(-9,0,0) +(-8,0,0) +(-7,0,0) +(-6,0,0) +(-12,6,0) +(-11,6,0) +(-10,6,0) +(-9,6,0) +(-8,6,0) +(-7,6,0) +(-6,6,0) +(-12,1,0) +(-12,2,0) +(-12,3,0) +(-12,4,0) +(-12,5,0) +(-6,1,0) +(-6,2,0) +(-6,3,0) +(-6,4,0) +(-6,5,0)