From: Jorge Arévalo Date: Wed, 16 Feb 2011 18:19:57 +0000 (+0000) Subject: - One raster core implementation of MapAlgebra (related ticket #588, needs X-Git-Tag: 2.0.0alpha1~1968 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=1fcfb1f7c0cee8d3512a480532a73ae7ecb0beed;p=postgis - One raster core implementation of MapAlgebra (related ticket #588, needs documentation). - RASTER_addBand code moved to core level. The new RASTER_addBand function calls the core one (rt_raster_generate_new_band). - Added regression tests for MapAlgebra. - Deleted lexer/parser at core level. Not used. - Fixed small bug in documentation: ST_SetBandNoDataValue returns a raster, not an integer. git-svn-id: http://svn.osgeo.org/postgis/trunk@6832 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index b50eb8295..bcfa38a52 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -2741,12 +2741,12 @@ WHERE rid = 2; - integer ST_SetBandNoDataValue + raster ST_SetBandNoDataValue raster rast double precision val - integer ST_SetBandNoDataValue + raster ST_SetBandNoDataValue raster rast integer band double precision val @@ -3713,4 +3713,4 @@ WHERE A.rid =2 ; - \ No newline at end of file + diff --git a/raster/rt_core/ma_lexer.l b/raster/rt_core/ma_lexer.l deleted file mode 100644 index 8d2e324bd..000000000 --- a/raster/rt_core/ma_lexer.l +++ /dev/null @@ -1,35 +0,0 @@ -%{ -#include -#include "ma_parser.tab.h" - -static YY_BUFFER_STATE ma_yy_buf_state; - - -/* -* Set up the lexer! -*/ -void ma_lexer_init(char *src) -{ - ma_yy_buf_state = yy_scan_string(src); -} - -/* -* Clean up the lexer! -*/ -void ma_lexer_close() -{ - yy_delete_buffer(ma_yy_buf_state); -} - -%} - -integer -?[0-9]+ -real -?(([0-9]+\.?)|([0-9]*\.?[0-9]+)([eE][-+]?[0-9]+)?) -symbol "+"|"-"|"*"|"/"|"^"|"("|")" -%% - -{integer} {yylval.integer = atoi(yytext); return TKN_INTEGER;} -{real} {yylval.real = atof(yytext); return TKN_REAL;} -{symbol} {return yytext[0];} -. ; -%% \ No newline at end of file diff --git a/raster/rt_core/ma_parser.y b/raster/rt_core/ma_parser.y deleted file mode 100644 index b7521ca06..000000000 --- a/raster/rt_core/ma_parser.y +++ /dev/null @@ -1,86 +0,0 @@ -%{ -#include -#include -#include -#include - -void yyerror(char *); -%} - -%union { - int integer; - double real; -} - -%token TKN_INTEGER -%token TKN_REAL -%type integer_exp -%type real_exp - -%left '-' '+' -%left '*' '/' -%left NEG -%right '^' - -%start calc - -%% - -calc: - | calc line - ; - -line: - integer_exp { printf("%d\n", $1); } - | real_exp { printf("%f\n", $1); } - ; - -integer_exp: - TKN_INTEGER { $$ = $1; } - | integer_exp '+' integer_exp { $$ = $1 + $3; } - | integer_exp '-' integer_exp { $$ = $1 - $3; } - | integer_exp '*' integer_exp { $$ = $1 * $3; } - | integer_exp '^' integer_exp { $$ = pow($1, $3); } - | '(' integer_exp ')' { $$ = $2; } - ; - -real_exp: - TKN_REAL { $$ = $1; } - | real_exp '+' real_exp { $$ = $1 + $3; } - | real_exp '-' real_exp { $$ = $1 - $3; } - | real_exp '*' real_exp { $$ = $1 * $3; } - | real_exp '/' real_exp { $$ = $1 / $3; } - | real_exp '^' real_exp { $$ = pow($1, $3); } - | '(' real_exp ')' { $$ = $2; } - ; - - -%% - -/* Just for testing purposes */ -int main(int argc, char **argv) -{ - char * testStr = (char*)calloc(10, sizeof(char)); - memcpy(testStr, argv[1], 10*sizeof(char)); - - ma_lexer_init(testStr); - - yyparse(); - - ma_lexer_close(); - - free(testStr); - - return 0; -} - -void yyerror(char *msg) -{ - fprintf(stderr, "Error: %s\n", msg); -} - -/* To compile: - bison -d ma_parser.y - flex ma_parser.l - gcc lex.yy.c ma_parser.tab.c -o ma_parser -lfl -lm -*/ \ No newline at end of file diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index e396f9bed..bc009a3e4 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -1382,8 +1382,6 @@ rt_raster_add_band(rt_context ctx, rt_raster raster, rt_band band, int index) { sizeof (rt_band)*(raster->numBands + 1) ); - RASTER_DEBUGF(4, "realloc returned %p", raster->bands); - if (!raster->bands) { ctx->err("Out of virtual memory " "reallocating band pointers"); @@ -1391,6 +1389,8 @@ rt_raster_add_band(rt_context ctx, rt_raster raster, rt_band band, int index) { return -1; } + RASTER_DEBUGF(4, "realloc returned %p", raster->bands); + for (i = 0; i <= raster->numBands; ++i) { if (i == index) { oldband = raster->bands[i]; @@ -1403,9 +1403,221 @@ rt_raster_add_band(rt_context ctx, rt_raster raster, rt_band band, int index) { } raster->numBands++; + + RASTER_DEBUGF(4, "now raster has %d bands", raster->numBands); + return index; } + +int32_t +rt_raster_generate_new_band(rt_context ctx, rt_raster raster, rt_pixtype pixtype, + double initialvalue, uint32_t hasnodata, double nodatavalue, int index) +{ + rt_band band = NULL; + int width = 0; + int height = 0; + int numval = 0; + int datasize = 0; + int oldnumbands = 0; + int numbands = 0; + void * mem = NULL; + int32_t checkvalint = 0; + uint32_t checkvaluint = 0; + double checkvaldouble = 0; + float checkvalfloat = 0; + int i; + + + assert(NULL != ctx); + assert(NULL != raster); + + /* Make sure index is in a valid range */ + oldnumbands = rt_raster_get_num_bands(ctx, raster); + if (index < 0) + index = 0; + else if (index > rt_raster_get_num_bands(ctx, raster) + 1) + index = rt_raster_get_num_bands(ctx, raster) + 1; + + /* Determine size of memory block to allocate and allocate it */ + width = rt_raster_get_width(ctx, raster); + height = rt_raster_get_height(ctx, raster); + numval = width * height; + datasize = rt_pixtype_size(ctx, pixtype) * numval; + + mem = (int *)ctx->alloc(datasize); + if (!mem) { + ctx->err("Could not allocate memory for band"); + return -1; + } + + if (fabs(initialvalue - 0.0) < FLT_EPSILON) + memset(mem, 0, datasize); + else { + switch (pixtype) + { + case PT_1BB: + { + uint8_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint8_t) initialvalue&0x01; + checkvalint = ptr[0]; + break; + } + case PT_2BUI: + { + uint8_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint8_t) initialvalue&0x03; + checkvalint = ptr[0]; + break; + } + case PT_4BUI: + { + uint8_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint8_t) initialvalue&0x0F; + checkvalint = ptr[0]; + break; + } + case PT_8BSI: + { + int8_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (int8_t) initialvalue; + checkvalint = ptr[0]; + break; + } + case PT_8BUI: + { + uint8_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint8_t) initialvalue; + checkvalint = ptr[0]; + break; + } + case PT_16BSI: + { + int16_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (int16_t) initialvalue; + checkvalint = ptr[0]; + break; + } + case PT_16BUI: + { + uint16_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint16_t) initialvalue; + checkvalint = ptr[0]; + break; + } + case PT_32BSI: + { + int32_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (int32_t) initialvalue; + checkvalint = ptr[0]; + break; + } + case PT_32BUI: + { + uint32_t *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (uint32_t) initialvalue; + checkvaluint = ptr[0]; + break; + } + case PT_32BF: + { + float *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = (float) initialvalue; + checkvalfloat = ptr[0]; + break; + } + case PT_64BF: + { + double *ptr = mem; + for (i = 0; i < numval; i++) + ptr[i] = initialvalue; + checkvaldouble = ptr[0]; + break; + } + default: + { + ctx->err("Unknown pixeltype %d", pixtype); + ctx->dealloc(mem); + return -1; + } + } + } + +#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION + /* Overflow checking */ + switch (pixtype) + { + case PT_1BB: + case PT_2BUI: + case PT_4BUI: + case PT_8BSI: + case PT_8BUI: + case PT_16BSI: + case PT_16BUI: + case PT_32BSI: + { + if (fabs(checkvalint - initialvalue) > FLT_EPSILON) + ctx->warn("Initial pixel value for %s band got truncated from %f to %d", + rt_pixtype_name(ctx, pixtype), + initialvalue, checkvalint); + break; + } + case PT_32BUI: + { + if (fabs(checkvaluint - initialvalue) > FLT_EPSILON) + ctx->warn("Initial pixel value for %s band got truncated from %f to %u", + rt_pixtype_name(ctx, pixtype), + initialvalue, checkvaluint); + break; + } + case PT_32BF: + { + /* For float, because the initial value is a double, + there is very often a difference between the desired value and the obtained one */ + if (fabs(checkvalfloat - initialvalue) > FLT_EPSILON) + ctx->warn("Initial pixel value for %s band got truncated from %f to %g", + rt_pixtype_name(ctx, pixtype), + initialvalue, checkvalfloat); + break; + } + case PT_64BF: + { + if (fabs(checkvaldouble - initialvalue) > FLT_EPSILON) + ctx->warn("Initial pixel value for %s band got truncated from %f to %g", + rt_pixtype_name(ctx, pixtype), + initialvalue, checkvaldouble); + break; + } + } +#endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */ + + band = rt_band_new_inline(ctx, width, height, pixtype, hasnodata, nodatavalue, mem); + if (! band) { + ctx->err("Could not add band to raster. Aborting"); + ctx->dealloc(mem); + return -1; + } + index = rt_raster_add_band(ctx, raster, band, index); + numbands = rt_raster_get_num_bands(ctx, raster); + if (numbands == oldnumbands || index == -1) { + ctx->err("Could not add band to raster. Aborting"); + rt_band_destroy(ctx, band); + } + + return index; +} + + void rt_raster_cell_to_geopoint(rt_context ctx, rt_raster raster, double x, double y, @@ -3319,9 +3531,9 @@ int rt_raster_has_no_band(rt_context ctx, rt_raster raster, int nband) { * @param ctx: context, for thread safety * @param raster1: raster to copy band to * @param raster2: raster to copy band from - * @param nband1: band index of destination raster - * @param nband2: band index of source raster - * @return The band index of the first raster where the new band is copied. + * @param nband1: band index of source raster + * @param nband2: band index of destination raster + * @return The band index of the second raster where the new band is copied. */ int32_t rt_raster_copy_band(rt_context ctx, rt_raster raster1, rt_raster raster2, int nband1, int nband2) @@ -3339,21 +3551,21 @@ int32_t rt_raster_copy_band(rt_context ctx, rt_raster raster1, } /* Check bands limits */ - if (nband1 < 1) - nband1 = 1; - else if (nband1 > raster1->numBands) - nband1 = raster1->numBands; + if (nband1 < 0) + nband1 = 0; + else if (nband1 >= raster1->numBands) + nband1 = raster1->numBands - 1; - if (nband2 < 1) - nband2 = 1; + if (nband2 < 0) + nband2 = 0; else if (nband2 > raster2->numBands) nband2 = raster2->numBands; - /* Get band from second raster */ - newband = raster2->bands[nband2]; + /* Get band from first raster */ + newband = rt_raster_get_band(ctx, raster1, nband1); - /* Add band to the first raster */ - return rt_raster_add_band(ctx, raster1, newband, nband1); + /* Add band to the second raster */ + return rt_raster_add_band(ctx, raster2, newband, nband2); } /** diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index e3b682d5c..8237c60ed 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -494,6 +494,24 @@ uint16_t rt_raster_get_height(rt_context ctx, rt_raster raster); */ int32_t rt_raster_add_band(rt_context ctx, rt_raster raster, rt_band band, int index); + + +/** + * Generate a new band data and add it to a raster. + * + * @param ctx : context, for thread safety + * @param raster : the raster to add a band to + * @param pixtype: the pixel type for the new band + * @param initialvalue: initial value for pixels + * @param hasnodata: indicates if the band has a nodata value + * @param nodatavalue: nodata value for the new band + * @param index: position to add the new band in the raster + * + * @return identifier (position) for the just-added raster, or -1 on error + */ +int32_t rt_raster_generate_new_band(rt_context ctx, rt_raster raster, rt_pixtype pixtype, + double initialvalue, uint32_t hasnodata, double nodatavalue, int index); + /** * Set scale in projection units * diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 6787a7547..05efbbf8e 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -38,6 +38,8 @@ #include #include #include +#include +#include #include /*#include "lwgeom_pg.h"*/ @@ -61,7 +63,9 @@ static rt_context get_rt_context(FunctionCallInfoData *fcinfo); static void *rt_pgalloc(size_t size); static void *rt_pgrealloc(void *mem, size_t size); static void rt_pgfree(void *mem); - +static char * replace(const char *str, const char *oldstr, const char *newstr, + int *count); +static char *strtoupper(char *str); /* Prototypes */ @@ -108,10 +112,79 @@ Datum RASTER_getBandPath(PG_FUNCTION_ARGS); Datum RASTER_getPixelValue(PG_FUNCTION_ARGS); Datum RASTER_setPixelValue(PG_FUNCTION_ARGS); Datum RASTER_addband(PG_FUNCTION_ARGS); +Datum RASTER_mapAlgebra(PG_FUNCTION_ARGS); Datum RASTER_isEmpty(PG_FUNCTION_ARGS); Datum RASTER_hasNoBand(PG_FUNCTION_ARGS); +/* Replace function taken from http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 */ +/* --------------------------------------------------------------------------- + Name : replace - Search & replace a substring by another one. + Creation : Thierry Husson, Sept 2010 + Parameters : + str : Big string where we search + oldstr : Substring we are looking for + newstr : Substring we want to replace with + count : Optional pointer to int (input / output value). NULL to ignore. + Input: Maximum replacements to be done. NULL or < 1 to do all. + Output: Number of replacements done or -1 if not enough memory. + Returns : Pointer to the new string or NULL if error. + Notes : + - Case sensitive - Otherwise, replace functions "strstr" by "strcasestr" + - Always allocates memory for the result. +--------------------------------------------------------------------------- */ +static char* +replace(const char *str, const char *oldstr, const char *newstr, int *count) +{ + const char *tmp = str; + char *result; + int found = 0; + int length, reslen; + int oldlen = strlen(oldstr); + int newlen = strlen(newstr); + int limit = (count != NULL && *count > 0) ? *count : -1; + + tmp = str; + while ((tmp = strstr(tmp, oldstr)) != NULL && found != limit) + found++, tmp += oldlen; + + length = strlen(str) + found * (newlen - oldlen); + if ( (result = (char *)palloc(length+1)) == NULL) { + fprintf(stderr, "Not enough memory\n"); + found = -1; + } else { + tmp = str; + limit = found; /* Countdown */ + reslen = 0; /* length of current result */ + /* Replace each old string found with new string */ + while ((limit-- > 0) && (tmp = strstr(tmp, oldstr)) != NULL) { + length = (tmp - str); /* Number of chars to keep intouched */ + strncpy(result + reslen, str, length); /* Original part keeped */ + strcpy(result + (reslen += length), newstr); /* Insert new string */ + reslen += newlen; + tmp += oldlen; + str = tmp; + } + strcpy(result + reslen, str); /* Copies last part and ending nul char */ + } + if (count != NULL) *count = found; + return result; +} + + +static char * +strtoupper(char * str) +{ + int j; + + for(j = 0; j < strlen(str); j++) + str[j] = toupper(str[j]); + + return str; + +} + + static void * rt_pgalloc(size_t size) { @@ -1281,7 +1354,7 @@ Datum RASTER_setBandNoDataValue(PG_FUNCTION_ARGS) bool forceChecking = FALSE; /* Check index is not NULL */ - if ((Pointer *)PG_GETARG_DATUM(1) == NULL) { + if (PG_ARGISNULL(1)) { /* Simply return NULL */ PG_RETURN_NULL(); } @@ -1295,7 +1368,7 @@ Datum RASTER_setBandNoDataValue(PG_FUNCTION_ARGS) assert(0 <= (index - 1)); /* Deserialize raster */ - if ((Pointer *)PG_GETARG_DATUM(0) == NULL) { + if (PG_ARGISNULL(0)) { /* Simply return NULL */ PG_RETURN_NULL(); } @@ -1319,7 +1392,7 @@ Datum RASTER_setBandNoDataValue(PG_FUNCTION_ARGS) forceChecking = PG_GETARG_BOOL(3); - if ((Pointer *)PG_GETARG_DATUM(2) == NULL) { + if (PG_ARGISNULL(2)) { /* Set the hasnodata flag to FALSE */ rt_band_set_hasnodata_flag(ctx, band, FALSE); @@ -1737,178 +1810,11 @@ Datum RASTER_addband(PG_FUNCTION_ARGS) index = oldnumbands + 1; } } - - /* Determine size of memory block to allocate and allocate*/ - width = rt_raster_get_width(ctx, raster); - height = rt_raster_get_height(ctx, raster); - numval = width * height; - datasize = rt_pixtype_size(ctx, pixtype) * numval; - - - mem = palloc((int) datasize); - if (!mem) - { - elog(ERROR, "Could not allocate memory for band"); - PG_RETURN_NULL(); - } - - /* Initialize block of memory to 0 or to the provided initial value */ - if (initialvalue == 0) - memset(mem, 0, datasize); - else - { - switch (pixtype) - { - case PT_1BB: - { - uint8_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint8_t) initialvalue&0x01; - checkvalint = ptr[0]; - break; - } - case PT_2BUI: - { - uint8_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint8_t) initialvalue&0x03; - checkvalint = ptr[0]; - break; - } - case PT_4BUI: - { - uint8_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint8_t) initialvalue&0x0F; - checkvalint = ptr[0]; - break; - } - case PT_8BSI: - { - int8_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (int8_t) initialvalue; - checkvalint = ptr[0]; - break; - } - case PT_8BUI: - { - uint8_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint8_t) initialvalue; - checkvalint = ptr[0]; - break; - } - case PT_16BSI: - { - int16_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (int16_t) initialvalue; - checkvalint = ptr[0]; - break; - } - case PT_16BUI: - { - uint16_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint16_t) initialvalue; - checkvalint = ptr[0]; - break; - } - case PT_32BSI: - { - int32_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (int32_t) initialvalue; - checkvalint = ptr[0]; - break; - } - case PT_32BUI: - { - uint32_t *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (uint32_t) initialvalue; - checkvaluint = ptr[0]; - break; - } - case PT_32BF: - { - float *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = (float) initialvalue; - checkvalfloat = ptr[0]; - break; - } - case PT_64BF: - { - double *ptr = mem; - for (i = 0; i < numval; i++) - ptr[i] = initialvalue; - checkvaldouble = ptr[0]; - break; - } - default: - { - elog(ERROR, "Unknown pixeltype %d", pixtype); - PG_RETURN_NULL(); - } - } - } -#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION - /* Overflow checking */ - switch (pixtype) - { - case PT_1BB: - case PT_2BUI: - case PT_4BUI: - case PT_8BSI: - case PT_8BUI: - case PT_16BSI: - case PT_16BUI: - case PT_32BSI: - { - if (fabs(checkvalint - initialvalue) > 0.0001) - elog(WARNING, "Initial pixel value for %s band got truncated from %f to %d", - rt_pixtype_name(ctx, pixtype), - initialvalue, checkvalint); - break; - } - case PT_32BUI: - { - if (fabs(checkvaluint - initialvalue) > 0.0001) - elog(WARNING, "Initial pixel value for %s band got truncated from %f to %u", - rt_pixtype_name(ctx, pixtype), - initialvalue, checkvaluint); - break; - } - case PT_32BF: - { - /* For float, because the initial value is a double, - there is very often a difference between the desired value and the obtained one */ - if (checkvalfloat != initialvalue) - elog(WARNING, "Initial pixel value for %s band got truncated from %f to %g", - rt_pixtype_name(ctx, pixtype), - initialvalue, checkvalfloat); - break; - } - case PT_64BF: - { - if (checkvaldouble != initialvalue) - elog(WARNING, "Initial pixel value for %s band got truncated from %f to %g", - rt_pixtype_name(ctx, pixtype), - initialvalue, checkvaldouble); - break; - } - } -#endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */ + + index = rt_raster_generate_new_band(ctx, raster, pixtype, initialvalue, + hasnodata, nodatavalue, index - 1); - band = rt_band_new_inline(ctx, width, height, pixtype, hasnodata, nodatavalue, mem); - if (! band) { - elog(ERROR, "Could not add band to raster. Returning NULL"); - PG_RETURN_NULL(); - } - index = rt_raster_add_band(ctx, raster, band, index - 1); numbands = rt_raster_get_num_bands(ctx, raster); if (numbands == oldnumbands || index == -1) { elog(ERROR, "Could not add band to raster. Returning NULL"); @@ -1975,6 +1881,445 @@ Datum RASTER_hasNoBand(PG_FUNCTION_ARGS) PG_RETURN_BOOL(rt_raster_has_no_band(ctx, raster, nBand)); } +PG_FUNCTION_INFO_V1(RASTER_mapAlgebra); +Datum RASTER_mapAlgebra(PG_FUNCTION_ARGS) +{ + rt_pgraster *pgraster = NULL; + rt_raster raster = NULL; + rt_raster newrast = NULL; + rt_context ctx = NULL; + rt_band band = NULL; + rt_band newband = NULL; + int x, y, nband, width, height; + double r; + double newnodatavalue, newinitialvalue, newval; + char *newexpr = NULL; + char *initexpr = NULL; + char *initndvexpr = NULL; + char *expression = NULL; + char *nodatavalueexpr = NULL; + rt_pixtype newpixeltype, pixeltype; + int skipcomputation = 0; + bool newhasnodatavalue = FALSE; + char strnewnodatavalue[50]; + char strnewval[50]; + int count = 0; + char *command = NULL; + int len = 0; + int ret = -1; + int proc; + TupleDesc tupdesc; + SPITupleTable * tuptable = NULL; + HeapTuple tuple; + + POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebra: STARTING..."); + + /* Check raster */ + if (PG_ARGISNULL(0)) { + elog(WARNING, "MapAlgebra: Raster is NULL. Returning NULL"); + PG_RETURN_NULL(); + } + + + /* Deserialize raster */ + pgraster = (rt_pgraster *)PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0)); + ctx = get_rt_context(fcinfo); + raster = rt_raster_deserialize(ctx, pgraster); + if ( ! raster ) + { + ereport(ERROR, + (errcode(ERRCODE_OUT_OF_MEMORY), + errmsg("MapAlgebra: Could not deserialize raster"))); + PG_RETURN_NULL(); + } + + POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebra: GETTING ARGUMENTS..."); + + /* Get the rest of the arguments */ + + if (PG_ARGISNULL(1)) + nband = 1; + else + nband = PG_GETARG_INT32(1); + + if (nband < 1) + nband = 1; + + POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebra: CREATING NEW RASTER..."); + + /** + * Create a new empty raster with having the same georeference as the + * provided raster + **/ + width = rt_raster_get_width(ctx, raster); + height = rt_raster_get_height(ctx, raster); + + newrast = rt_raster_new(ctx, + rt_raster_get_width(ctx, raster), + rt_raster_get_height(ctx, raster)); + + if ( ! newrast ) { + elog(ERROR, "MapAlgebra: Could not create a new raster"); + PG_RETURN_NULL(); + } + + rt_raster_set_scale(ctx, newrast, + rt_raster_get_x_scale(ctx, raster), + rt_raster_get_y_scale(ctx, raster)); + + rt_raster_set_offsets(ctx, newrast, + rt_raster_get_x_offset(ctx, raster), + rt_raster_get_y_offset(ctx, raster)); + + rt_raster_set_skews(ctx, newrast, + rt_raster_get_x_skew(ctx, raster), + rt_raster_get_y_skew(ctx, raster)); + + rt_raster_set_srid(ctx, newrast, rt_raster_get_srid(ctx, raster)); + + + + /** + * If this new raster is empty (width = 0 OR height = 0) then there is + * nothing to compute and we return it right now + **/ + if (rt_raster_is_empty(ctx, newrast)) { + elog(WARNING, "MapAlgebra: Raster is empty. Returning an empty raster"); + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + PG_RETURN_POINTER(pgraster); + } + + POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebra: GETTING RASTER BAND %d...", nband); + + + /** + * Check if the raster has the required band. Otherwise, return a raster + * without band + **/ + if (rt_raster_has_no_band(ctx, raster, nband)) { + elog(WARNING, "Mapalgebra: Raster do not have the required band. " + "Returning a raster without a band"); + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + PG_RETURN_POINTER(pgraster); + } + + /* Get the raster band */ + band = rt_raster_get_band(ctx, raster, nband - 1); + if ( ! band ) { + elog(ERROR, "MapAlgebra: Could not get band %d for raster", nband); + PG_RETURN_NULL(); + } + + /* Check for nodata value */ + if (rt_band_get_hasnodata_flag(ctx, band)) { + newnodatavalue = rt_band_get_nodata(ctx, band); + } + + else { + newnodatavalue = rt_band_get_min_value(ctx, band); + } + + POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebra: NODATA VALUE FOR BAND: %f...", + newnodatavalue); + /** + * We set the initial value of the future band to nodata value. If nodata + * value is null, then the raster will be initialized to + * rt_band_get_min_value but all the values should be recomputed anyway + **/ + newinitialvalue = newnodatavalue; + + POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebra: SETTING PIXELTYPE..."); + + /** + * Set the new pixeltype + **/ + if (PG_ARGISNULL(4)) { + newpixeltype = rt_band_get_pixtype(ctx, band); + } + + else { + newpixeltype = rt_pixtype_index_from_name(ctx, + text_to_cstring(PG_GETARG_TEXT_P(4))); + if (newpixeltype == PT_END) + newpixeltype = rt_band_get_pixtype(ctx, band); + } + + if (newpixeltype == PT_END) { + elog(ERROR, "MapAlgebra: Invalid pixeltype. Aborting"); + PG_RETURN_NULL(); + } + + /* Connect with SPI manager */ + SPI_connect(); + + /* Construct expression for raster values */ + + if (!PG_ARGISNULL(2)) { + expression = text_to_cstring(PG_GETARG_TEXT_P(2)); + len = strlen("SELECT ") + strlen(expression); + initexpr = (char *)palloc(len + 1); + + strncpy(initexpr, "SELECT ", strlen("SELECT ")); + strncpy(initexpr + strlen("SELECT "), strtoupper(expression), + strlen(expression)); + initexpr[len] = '\0'; + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: Expression is %s", initexpr); + + } + + /** + * Optimization: If a nodatavalueexpr is provided, recompute the initial + * value. Then, we can initialize the raster with this value and skip the + * computation of nodata values one by one in the main computing loop + **/ + if (!PG_ARGISNULL(3)) { + nodatavalueexpr = text_to_cstring(PG_GETARG_TEXT_P(3)); + len = strlen("SELECT ") + strlen(nodatavalueexpr); + initndvexpr = (char *)palloc(len + 1); + strncpy(initndvexpr, "SELECT ", strlen("SELECT ")); + strncpy(initndvexpr + strlen("SELECT "), strtoupper(nodatavalueexpr), + strlen(nodatavalueexpr)); + initndvexpr[len] = '\0'; + sprintf(strnewnodatavalue, "%f", newnodatavalue); + + newexpr = replace(initndvexpr, "RAST", strnewnodatavalue, &count); + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: initndvexpr = %s", initndvexpr); + + /** + * Execute the expression for nodata value and store the result as new + * initial value + **/ + ret = SPI_execute(newexpr, FALSE, 0); + + if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) { + elog(ERROR, "MapAlgebra: invalid construction for NODATA " + "expression. Aborting"); + SPI_finish(); + PG_RETURN_NULL(); + } + + tupdesc = SPI_tuptable->tupdesc; + tuptable = SPI_tuptable; + + tuple = tuptable->vals[0]; + newinitialvalue = atof(SPI_getvalue(tuple, tupdesc, 1)); + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: new initial value = %f", + newinitialvalue); + + } + + + /** + * Optimization: If the raster is only filled with nodata values return + * right now a raster filled with the nodatavalueexpr + * TODO: Call rt_band_check_isnodata instead? + **/ + if (rt_band_get_isnodata_flag(ctx, band)) { + + POSTGIS_RT_DEBUG(3, "MapAlgebra: Band is a nodata band, returning " + "a raster filled with nodata"); + + ret = rt_raster_generate_new_band(ctx, newrast, newpixeltype, + newinitialvalue, TRUE, newnodatavalue, 0); + + /* Serialize created raster */ + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + + /* Disconnect function from SPI manager */ + SPI_finish(); + + PG_RETURN_POINTER(pgraster); + } + + + /** + * Optimization: If expression resume to 'RAST' and nodatavalueexpr is NULL + * or also equal to 'RAST', we can just return the band from the original + * raster + **/ + if (initexpr != NULL && !strcmp(initexpr, "SELECT RAST") && + (nodatavalueexpr == NULL || !strcmp(initndvexpr, "SELECT RAST"))) { + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: Expression resumes to RAST. Returning " + "raster with band %d from original raster", nband); + + POSTGIS_RT_DEBUGF(4, "MapAlgebra: New raster has %d bands", + rt_raster_get_num_bands(ctx, newrast)); + + rt_raster_copy_band(ctx, raster, newrast, nband - 1, 0); + + POSTGIS_RT_DEBUGF(4, "MapAlgebra: New raster now has %d bands", + rt_raster_get_num_bands(ctx, newrast)); + + /* Serialize created raster */ + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + + /* Disconnect function from SPI manager */ + SPI_finish(); + + PG_RETURN_POINTER(pgraster); + } + + /** + * Optimization: If expression resume to a constant (it does not contain + * RAST) + **/ + if (initexpr != NULL && strstr(initexpr, "RAST") == NULL) { + /* Execute the expresion into newval */ + ret = SPI_execute(initexpr, FALSE, 0); + + if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) { + elog(ERROR, "MapAlgebra: invalid construction for expression. Aborting"); + SPI_finish(); + PG_RETURN_NULL(); + } + + tupdesc = SPI_tuptable->tupdesc; + tuptable = SPI_tuptable; + + tuple = tuptable->vals[0]; + newval = atof(SPI_getvalue(tuple, tupdesc, 1)); + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: new raster value = %f", + newval); + + skipcomputation = 1; + + /** + * Compute the new value, set it and we will return after creating the + * new raster + **/ + if (nodatavalueexpr == NULL) { + newinitialvalue = newval; + skipcomputation = 2; + } + + /* Return the new raster as it will be before computing pixel by pixel */ + else if (fabs(newval - newinitialvalue) > FLT_EPSILON) { + skipcomputation = 2; + } + } + + /** + * Create the raster receiving all the computed values. Initialize it to the + * new initial value + **/ + ret = rt_raster_generate_new_band(ctx, newrast, newpixeltype, + newinitialvalue, TRUE, newnodatavalue, 0); + + /** + * Optimization: If expression is NULL, or all the pixels could be set in + * one step, return the initialized raster now + **/ + if (expression == NULL || skipcomputation == 2) { + /* Serialize created raster */ + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + + /* Disconnect function from SPI manager */ + SPI_finish(); + + PG_RETURN_POINTER(pgraster); + } + + /* Get the new raster band */ + newband = rt_raster_get_band(ctx, newrast, 0); + if ( ! newband ) { + elog(WARNING, "MapAlgebra: Could not modify band for new raster. " + "Returning new raster with the original band"); + + /* Serialize created raster */ + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + + /* Disconnect function from SPI manager */ + SPI_finish(); + + PG_RETURN_POINTER(pgraster); + } + + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: Main computing loop (%d x %d)", + width, height); + + for (x = 0; x < width; x++) { + for(y = 0; y < height; y++) { + ret = rt_band_get_pixel(ctx, band, x, y, &r); + + /** + * We compute a value only for the withdata value pixel since the + * nodata value has already been set by the first optimization + **/ + if (ret != -1 && fabs(r - newnodatavalue) > FLT_EPSILON) { + if (skipcomputation == 0) { + sprintf(strnewval, "%f", r); + + if (initexpr != NULL) { + newexpr = replace(initexpr, "RAST", strnewval, &count); + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: (%dx%d), r = %s, newexpr = %s", + x, y, strnewval, newexpr); + + ret = SPI_execute(newexpr, FALSE, 0); + + if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) { + elog(ERROR, "MapAlgebra: invalid construction for expression. Aborting"); + SPI_finish(); + PG_RETURN_NULL(); + } + + tupdesc = SPI_tuptable->tupdesc; + tuptable = SPI_tuptable; + + tuple = tuptable->vals[0]; + newval = atof(SPI_getvalue(tuple, tupdesc, 1)); + } + + else + newval = newinitialvalue; + + POSTGIS_RT_DEBUGF(3, "MapAlgebra: new value = %f", newval); + } + + + rt_band_set_pixel(ctx, newband, x, y, newval); + } + + } + } + + /* The newrast band has been modified */ + + /* Serialize created raster */ + pgraster = rt_raster_serialize(ctx, newrast); + if (!pgraster) PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + + /* Disconnect function from SPI manager */ + SPI_finish(); + + PG_RETURN_POINTER(pgraster); +} /* ---------------------------------------------------------------- */ /* Memory allocation / error reporting hooks */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 8de8e75c0..383c6859c 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -209,8 +209,54 @@ CREATE OR REPLACE FUNCTION st_addband(rast raster, index int, pixeltype text, in RETURNS raster AS 'select st_addband($1, $2, $3, $4, NULL)' LANGUAGE 'SQL' IMMUTABLE; + + +----------------------------------------------------------------------- +-- MapAlgebra +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, band integer, + expression text, nodatavalueexpr text, pixeltype text) + RETURNS raster + AS 'MODULE_PATHNAME', 'RASTER_mapAlgebra' + LANGUAGE 'C' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, band integer, + expression text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, $2, $3, NULL, NULL) $$ + LANGUAGE SQL; + +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, expression text, + pixeltype text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, 1, $2, NULL, $3) $$ + LANGUAGE SQL; + +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, expression text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, 1, $2, NULL, NULL) $$ + LANGUAGE SQL; + +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, band integer, + expression text, nodatavalueexpr text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, $2, $3, $4, NULL) $$ + LANGUAGE SQL; + + +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, expression text, + nodatavalueexpr text, pixeltype text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, 1, $2, $3, $4) $$ + LANGUAGE SQL; - +CREATE OR REPLACE FUNCTION st_mapalgebra(rast raster, expression text, + nodatavalueexpr text) + RETURNS raster + AS $$ SELECT st_mapalgebra($1, 1, $2, $3, NULL) $$ + LANGUAGE SQL; + + ----------------------------------------------------------------------- -- Get information about the raster ----------------------------------------------------------------------- diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 69a6b4f6c..afabbe12e 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -53,6 +53,8 @@ TEST_PIXEL = \ TEST_UTILITY = \ create_rt_utility_test.sql \ rt_utility.sql \ + create_rt_mapalgebra_test.sql \ + rt_mapalgebra.sql \ $(NULL) TEST_GIST = \ diff --git a/raster/test/regress/create_rt_mapalgebra_test.sql b/raster/test/regress/create_rt_mapalgebra_test.sql new file mode 100644 index 000000000..b1d0d8afd --- /dev/null +++ b/raster/test/regress/create_rt_mapalgebra_test.sql @@ -0,0 +1,11 @@ +CREATE OR REPLACE FUNCTION ST_TestRaster(ulx float8, uly float8, val float8) + RETURNS raster AS + $$ + DECLARE + BEGIN + RETURN ST_AddBand(ST_MakeEmptyRaster(10, 10, ulx, uly, 1, 1, 0, 0, -1), '32BF', val, -1); + END; + $$ + LANGUAGE 'plpgsql'; + + diff --git a/raster/test/regress/create_rt_mapalgebra_test_expected b/raster/test/regress/create_rt_mapalgebra_test_expected new file mode 100644 index 000000000..e69de29bb diff --git a/raster/test/regress/rt_mapalgebra.sql b/raster/test/regress/rt_mapalgebra.sql new file mode 100644 index 000000000..d39b3bbc8 --- /dev/null +++ b/raster/test/regress/rt_mapalgebra.sql @@ -0,0 +1,23 @@ +-- Test NULL raster +SELECT ST_MapAlgebra(NULL, 1, 'rast + 20', '2') IS NULL FROM ST_TestRaster(0, 0, -1) rast; + +-- Test empty raster +SELECT ST_IsEmpty(ST_MapAlgebra(ST_MakeEmptyRaster(0, 10, 0, 0, 1, 1, 1, 1, -1), 1, 'rast + 20', '2')); + +-- Test hasnoband raster +SELECT ST_HasNoBand(ST_MapAlgebra(ST_MakeEmptyRaster(10, 10, 0, 0, 1, 1, 1, 1, -1), 1, 'rast + 20', '2')); + +-- Test hasnodata value +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(ST_SetBandNoDataValue(rast, NULL), 1, 'rast + 20', '2'), 1, 1) FROM ST_TestRaster(0, 0, -1) rast; + +-- Test nodata value +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, 'rast + 20', 'rast + 2'), 1, 1) FROM ST_TestRaster(0, 0, -1) rast; + +-- Test 'rast' expression +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, 'rast', 'rast'), 1, 1) FROM ST_TestRaster(0, 0, -1) rast; +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(ST_SetBandNoDataValue(rast, NULL), 1, 'rast', NULL), 1, 1) FROM ST_TestRaster(0, 0, -1) rast; + +-- Test pixeltype +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, 'rast + 20', 'rast + 2', '4BUI'), 1, 1) FROM ST_TestRaster(0, 0, 100) rast; +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, 'rast + 20', 'rast + 2', '4BUId'), 1, 1) FROM ST_TestRaster(0, 0, 100) rast; +SELECT ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, 'rast + 20', 'rast + 2', '2BUI'), 1, 1) FROM ST_TestRaster(0, 0, 101) rast; diff --git a/raster/test/regress/rt_mapalgebra_expected b/raster/test/regress/rt_mapalgebra_expected new file mode 100644 index 000000000..f8c810dbc --- /dev/null +++ b/raster/test/regress/rt_mapalgebra_expected @@ -0,0 +1,213 @@ +WARNING: MapAlgebra: Raster is NULL. Returning NULL +t +WARNING: MapAlgebra: Raster is empty. Returning an empty raster +t +WARNING: Mapalgebra: Raster do not have the required band. Returning a raster without a band +t +|19 +|1 +| +|-1 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +WARNING: Pixel value for 4BUI band got truncated from 120 to 8 +100|8 +100|120 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +WARNING: Pixel value for 2BUI band got truncated from 121 to 1 +101|1