From c4d679a6d6cfc7e6c004496f258b3e832ad5815a Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Fri, 2 Mar 2012 01:15:35 +0000 Subject: [PATCH] Add handling of when clamped pixel values is equal to the band's clamped NODATA value. Ticket is #1139. git-svn-id: http://svn.osgeo.org/postgis/trunk@9380 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/rt_core/rt_api.c | 422 ++++++++++++++++++++++++++----------- raster/rt_core/rt_api.h | 27 ++- raster/rt_pg/rt_pg.c | 2 +- raster/test/core/testapi.c | 61 ------ 4 files changed, 322 insertions(+), 190 deletions(-) diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index b1166768a..f06c8124f 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -1804,145 +1804,145 @@ rt_band_set_pixel_line( * 1 on truncation/clamping/converting. */ int -rt_band_set_pixel(rt_band band, uint16_t x, uint16_t y, - double val) { - rt_pixtype pixtype = PT_END; - unsigned char* data = NULL; - uint32_t offset = 0; - - int32_t checkvalint = 0; - uint32_t checkvaluint = 0; - float checkvalfloat = 0; - double checkvaldouble = 0; - - double checkval = 0; - +rt_band_set_pixel( + rt_band band, + uint16_t x, uint16_t y, + double val +) { + rt_pixtype pixtype = PT_END; + unsigned char* data = NULL; + uint32_t offset = 0; + int rtn = 0; + int32_t checkvalint = 0; + uint32_t checkvaluint = 0; + float checkvalfloat = 0; + double checkvaldouble = 0; + double checkval = 0; - assert(NULL != band); + assert(NULL != band); - if (band->offline) { - rterror("rt_band_set_pixel not implemented yet for OFFDB bands"); - return -1; - } + if (band->offline) { + rterror("rt_band_set_pixel not implemented yet for OFFDB bands"); + return -1; + } - pixtype = band->pixtype; + pixtype = band->pixtype; - if (x >= band->width || y >= band->height) { - rterror("rt_band_set_pixel: Coordinates out of range"); - return -1; - } + if (x >= band->width || y >= band->height) { + rterror("rt_band_set_pixel: Coordinates out of range"); + return -1; + } - data = rt_band_get_data(band); - offset = x + (y * band->width); + /* check that clamped value isn't clamped NODATA */ + if (band->hasnodata && pixtype != PT_64BF) { + double newval; + if (rt_band_corrected_clamped_value(band, val, &newval) == 1) { +#if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0 + rtwarn("Value for pixel %d x %d has been corrected as clamped value becomes NODATA", x, y); +#endif + val = newval; + rtn = 1; + } + } - switch (pixtype) { - case PT_1BB: - { - data[offset] = rt_util_clamp_to_1BB(val); - checkvalint = data[offset]; - break; - } - case PT_2BUI: - { - data[offset] = rt_util_clamp_to_2BUI(val); - checkvalint = data[offset]; - break; - } - case PT_4BUI: - { - data[offset] = rt_util_clamp_to_4BUI(val); - checkvalint = data[offset]; - break; - } - case PT_8BSI: - { - data[offset] = rt_util_clamp_to_8BSI(val); - checkvalint = (int8_t) data[offset]; - break; - } - case PT_8BUI: - { - data[offset] = rt_util_clamp_to_8BUI(val); - checkvalint = data[offset]; - break; - } - case PT_16BSI: - { - int16_t *ptr = (int16_t*) data; /* we assume correct alignment */ - ptr[offset] = rt_util_clamp_to_16BSI(val); - checkvalint = (int16_t) ptr[offset]; - break; - } - case PT_16BUI: - { - uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */ - ptr[offset] = rt_util_clamp_to_16BUI(val); - checkvalint = ptr[offset]; - break; - } - case PT_32BSI: - { - int32_t *ptr = (int32_t*) data; /* we assume correct alignment */ - ptr[offset] = rt_util_clamp_to_32BSI(val); - checkvalint = (int32_t) ptr[offset]; - break; - } - case PT_32BUI: - { - uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */ - ptr[offset] = rt_util_clamp_to_32BUI(val); - checkvaluint = ptr[offset]; - break; - } - case PT_32BF: - { - float *ptr = (float*) data; /* we assume correct alignment */ - ptr[offset] = rt_util_clamp_to_32F(val); - checkvalfloat = ptr[offset]; - break; - } - case PT_64BF: - { - double *ptr = (double*) data; /* we assume correct alignment */ - ptr[offset] = val; - checkvaldouble = ptr[offset]; - break; - } - default: - { - rterror("rt_band_set_pixel: Unknown pixeltype %d", pixtype); - return -1; - } - } + data = rt_band_get_data(band); + offset = x + (y * band->width); - /* If the stored value is different from no data, reset the isnodata flag */ - if (FLT_NEQ(checkval, band->nodataval)) { - band->isnodata = FALSE; - } + switch (pixtype) { + case PT_1BB: { + data[offset] = rt_util_clamp_to_1BB(val); + checkvalint = data[offset]; + break; + } + case PT_2BUI: { + data[offset] = rt_util_clamp_to_2BUI(val); + checkvalint = data[offset]; + break; + } + case PT_4BUI: { + data[offset] = rt_util_clamp_to_4BUI(val); + checkvalint = data[offset]; + break; + } + case PT_8BSI: { + data[offset] = rt_util_clamp_to_8BSI(val); + checkvalint = (int8_t) data[offset]; + break; + } + case PT_8BUI: { + data[offset] = rt_util_clamp_to_8BUI(val); + checkvalint = data[offset]; + break; + } + case PT_16BSI: { + int16_t *ptr = (int16_t*) data; /* we assume correct alignment */ + ptr[offset] = rt_util_clamp_to_16BSI(val); + checkvalint = (int16_t) ptr[offset]; + break; + } + case PT_16BUI: { + uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */ + ptr[offset] = rt_util_clamp_to_16BUI(val); + checkvalint = ptr[offset]; + break; + } + case PT_32BSI: { + int32_t *ptr = (int32_t*) data; /* we assume correct alignment */ + ptr[offset] = rt_util_clamp_to_32BSI(val); + checkvalint = (int32_t) ptr[offset]; + break; + } + case PT_32BUI: { + uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */ + ptr[offset] = rt_util_clamp_to_32BUI(val); + checkvaluint = ptr[offset]; + break; + } + case PT_32BF: { + float *ptr = (float*) data; /* we assume correct alignment */ + ptr[offset] = rt_util_clamp_to_32F(val); + checkvalfloat = ptr[offset]; + break; + } + case PT_64BF: { + double *ptr = (double*) data; /* we assume correct alignment */ + ptr[offset] = val; + checkvaldouble = ptr[offset]; + break; + } + default: { + rterror("rt_band_set_pixel: Unknown pixeltype %d", pixtype); + return -1; + } + } - /* - * If the pixel was a nodata value, now the band may be NODATA band) - * TODO: NO, THAT'S TOO SLOW!!! - */ + /* If the stored value is different from no data, reset the isnodata flag */ + if (FLT_NEQ(checkval, band->nodataval)) { + band->isnodata = FALSE; + } - /* - else { - rt_band_check_is_nodata(band); - } - */ + /* + * If the pixel was a nodata value, now the band may be NODATA band) + * TODO: NO, THAT'S TOO SLOW!!! + */ + /* + else { + rt_band_check_is_nodata(band); + } + */ - /* Overflow checking */ - if (rt_util_dbl_trunc_warning( - val, - checkvalint, checkvaluint, - checkvalfloat, checkvaldouble, - pixtype - )) { - return 1; - } + /* Overflow checking */ + if (rt_util_dbl_trunc_warning( + val, + checkvalint, checkvaluint, + checkvalfloat, checkvaldouble, + pixtype + )) { + return 1; + } - return 0; + return rtn; } /** @@ -2139,6 +2139,174 @@ rt_band_check_is_nodata(rt_band band) return TRUE; } +/** + * Compare clamped value to band's clamped NODATA value. If unclamped + * value is exactly unclamped NODATA value, function returns -1. + * + * @param band: the band whose NODATA value will be used for comparison + * @param val: the value to compare to the NODATA value + * + * @return 1 if clamped value is clamped NODATA + * 0 if clamped value is NOT clamped NODATA + * -1 otherwise + */ +int +rt_band_clamped_value_is_nodata(rt_band band, double val) { + + assert(NULL != band); + + /* no NODATA, so no need to test */ + if (!band->hasnodata) + return -1; + + /* value is exactly NODATA */ + if (FLT_EQ(val, band->nodataval)) + return -1; + + switch (band->pixtype) { + case PT_1BB: + if (rt_util_clamp_to_1BB(val) == rt_util_clamp_to_1BB(band->nodataval)) + return 1; + break; + case PT_2BUI: + if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(band->nodataval)) + return 1; + break; + case PT_4BUI: + if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(band->nodataval)) + return 1; + break; + case PT_8BSI: + if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(band->nodataval)) + return 1; + break; + case PT_8BUI: + if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(band->nodataval)) + return 1; + break; + case PT_16BSI: + if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(band->nodataval)) + return 1; + break; + case PT_16BUI: + if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(band->nodataval)) + return 1; + break; + case PT_32BSI: + if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(band->nodataval)) + return 1; + break; + case PT_32BUI: + if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(band->nodataval)) + return 1; + break; + case PT_32BF: + if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(band->nodataval))) + return 1; + break; + case PT_64BF: + break; + default: + rterror("rt_band_clamped_value_is_nodata: Unknown pixeltype %d", band->pixtype); + return -1; + } + + return 0; +} + +/** + * Correct value when clamped value is clamped NODATA value. Correction + * does NOT occur if unclamped value is exactly unclamped NODATA value. + * + * @param band: the band whose NODATA value will be used for comparison + * @param val: the value to compare to the NODATA value and correct + * @param newval: pointer to corrected value + * + * @return 0 on error, 1 if corrected, -1 otherwise + */ +int +rt_band_corrected_clamped_value(rt_band band, double val, double *newval) { + double minval = 0.; + + assert(NULL != band); + + /* check that value needs correcting */ + if (rt_band_clamped_value_is_nodata(band, val) != 1) { + *newval = val; + return -1; + } + + minval = rt_pixtype_get_min_value(band->pixtype); + *newval = val; + + switch (band->pixtype) { + case PT_1BB: + *newval = !band->nodataval; + break; + case PT_2BUI: + if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_4BUI: + if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_8BSI: + if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_8BUI: + if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_16BSI: + if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_16BUI: + if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_32BSI: + if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_32BUI: + if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(minval)) + (*newval)++; + else + (*newval)--; + break; + case PT_32BF: + if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(minval))) + *newval += FLT_EPSILON; + else + *newval -= FLT_EPSILON; + break; + case PT_64BF: + break; + default: + rterror("rt_band_alternative_clamped_value: Unknown pixeltype %d", band->pixtype); + return 0; + } + + return 1; +} + /** * Compute summary statistics for a band * diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 1e5c30cc4..c33641cb8 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -92,8 +92,9 @@ #include /* For size_t, srand and rand */ #include /* For C99 int types */ -#include /* for FLT_EPSILON and float type limits */ +#include /* for FLT_EPSILON, DBL_EPSILON and float type limits */ #include /* for integer type limits */ +#include #include "lwgeom_geos.h" #include "liblwgeom.h" @@ -551,6 +552,30 @@ int rt_band_is_nodata(rt_band band); */ int rt_band_check_is_nodata(rt_band band); +/** + * Compare clamped value to band's clamped NODATA value + * + * @param band: the band whose NODATA value will be used for comparison + * @param val: the value to compare to the NODATA value + * + * @return 1 if clamped value is clamped NODATA + * 0 if clamped value is NOT clamped NODATA + * -1 otherwise + */ +int rt_band_clamped_value_is_nodata(rt_band band, double val); + +/** + * Correct value when clamped value is clamped NODATA value. Correction + * does NOT occur if unclamped value is exactly unclamped NODATA value. + * + * @param band: the band whose NODATA value will be used for comparison + * @param val: the value to compare to the NODATA value and correct + * @param newval: pointer to corrected value + * + * @return 0 on error, 1 if corrected, -1 otherwise + */ +int +rt_band_corrected_clamped_value(rt_band band, double val, double *newval); /** * Compute summary statistics for a band diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index c5691c4cb..7acb46d40 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -4,7 +4,7 @@ * WKTRaster - Raster Types for PostGIS * http://www.postgis.org/support/wiki/index.php?WKTRasterHomePage * - * Copyright (C) 2011 Regents of the University of California + * Copyright (C) 2011-2012 Regents of the University of California * * Copyright (C) 2010-2011 Jorge Arevalo * Copyright (C) 2010-2011 David Zwarg diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 28e4ca808..5dda62c7d 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1,8 +1,6 @@ #include #include #include -#include /* for fabs */ -#include /* for FLT_EPSILON and DBL_EPSILON */ #include "rt_api.h" #include "check.h" @@ -59,36 +57,11 @@ fillRasterToPolygonize(int hasnodata, double nodatavalue) { /* Create raster */ - /* First test raster */ - /* - uint16_t width = 2; - uint16_t height = 2; - */ - - /* Second test raster */ uint16_t width = 9; uint16_t height = 9; - /* Third test raster */ - /* - uint16_t width = 5; - uint16_t height = 5; - */ - rt_raster raster = rt_raster_new(width, height); - /* Fill raster. Option 1: simple raster */ - /* - rt_band band = addBand(raster, PT_32BSI, 0, 0); - - rt_band_set_pixel(band, 0, 0, 1); - rt_band_set_pixel(band, 0, 1, 1); - rt_band_set_pixel(band, 1, 0, 1); - rt_band_set_pixel(band, 1, 1, 1); - */ - - - /* Fill raster. Option 2: 9x9, 1 band */ rt_band band = addBand(raster, PT_32BF, hasnodata, nodatavalue); { @@ -98,7 +71,6 @@ fillRasterToPolygonize(int hasnodata, double nodatavalue) rt_band_set_pixel(band, x, y, 0.0); } - /**/ rt_band_set_pixel(band, 3, 1, 1.8); rt_band_set_pixel(band, 4, 1, 1.8); rt_band_set_pixel(band, 5, 1, 2.8); @@ -128,39 +100,6 @@ fillRasterToPolygonize(int hasnodata, double nodatavalue) rt_band_set_pixel(band, 4, 7, 1.8); rt_band_set_pixel(band, 5, 7, 2.8); - - - /* Fill raster. Option 3: 5x5, 1 band */ - /* - rt_band band = addBand(raster, PT_8BUI, 1, 255); - - rt_band_set_pixel(band, 0, 0, 253); - rt_band_set_pixel(band, 1, 0, 254); - rt_band_set_pixel(band, 2, 0, 253); - rt_band_set_pixel(band, 3, 0, 254); - rt_band_set_pixel(band, 4, 0, 254); - rt_band_set_pixel(band, 0, 1, 253); - rt_band_set_pixel(band, 1, 1, 254); - rt_band_set_pixel(band, 2, 1, 254); - rt_band_set_pixel(band, 3, 1, 253); - rt_band_set_pixel(band, 4, 1, 249); - rt_band_set_pixel(band, 0, 2, 250); - rt_band_set_pixel(band, 1, 2, 254); - rt_band_set_pixel(band, 2, 2, 254); - rt_band_set_pixel(band, 3, 2, 252); - rt_band_set_pixel(band, 4, 2, 249); - rt_band_set_pixel(band, 0, 3, 251); - rt_band_set_pixel(band, 1, 3, 253); - rt_band_set_pixel(band, 2, 3, 254); - rt_band_set_pixel(band, 3, 3, 254); - rt_band_set_pixel(band, 4, 3, 253); - rt_band_set_pixel(band, 0, 4, 252); - rt_band_set_pixel(band, 1, 4, 250); - rt_band_set_pixel(band, 2, 4, 254); - rt_band_set_pixel(band, 3, 4, 254); - rt_band_set_pixel(band, 4, 4, 254); - */ - return raster; } -- 2.40.0