Add handling of when clamped pixel values is equal to the band's clamped NODATA value...

git-svn-id: http://svn.osgeo.org/postgis/trunk@9380 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index b1166768a0be5505211f59e674d2472676c3d3e0..f06c8124fcae8963b085dae9ca1c448e98be7404 100644 (file)
--- 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 1e5c30cc4825496174c51744da4f65a2f1d2e70b..c33641cb807ac5d6d90837e396da7a157feca0b9 100644 (file)
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -92,8 +92,9 @@
 
 #include <stdlib.h> /* For size_t, srand and rand */
 #include <stdint.h> /* For C99 int types */
-#include <float.h> /* for FLT_EPSILON and float type limits */
+#include <float.h> /* for FLT_EPSILON, DBL_EPSILON and float type limits */
 #include <limits.h> /* for integer type limits */
+#include <math.h>
 
 #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 c5691c4cbf468a2211127ce7d3fffd8c8957903e..7acb46d408efc3f7a0ccaf9bae6547a11c1e9858 100644 (file)
--- 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
  *   <bkpark@ucdavis.edu>
  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>

diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c
index 28e4ca808d9eec66cf8737a487994a27d773d2ff..5dda62c7d2f962171ed15010c05b847b2073a64d 100644 (file)
--- a/raster/test/core/testapi.c
+++ b/raster/test/core/testapi.c
@@ -1,8 +1,6 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <assert.h>
-#include <math.h> /* for fabs */
-#include <float.h> /* 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;
 }