New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PROCESSING SCALE doesn't scale colors #4927

Closed
HuidaeCho opened this Issue May 17, 2014 · 3 comments

Comments

Projects
None yet
3 participants
@HuidaeCho
Copy link

HuidaeCho commented May 17, 2014

PROCESSING "SCALE=*" scales pixel values, but doesn't scale the corresponding colors. For example, originally

pixel   color
0       10,10,10
100     0,0,100
1000    0,100,0
2500    100,0,0

If I use SCALE=AUTO, the pixel values are scaled such that they are from 0 to 255, but the colors are not scaled:

pixel   color
0       10,10,10
10      0,0,0 # it has to be 0,0,100 because pixel value 100 was scaled to 10.
100     0,0,100 # it has to be 0,100,0.
250     0,0,0 # it has to be 100,0,0.

Because of the above issue, the map is rendered differently from the original map. Expected results are:

pixel   color
0       10,10,10
10      0,0,100
100     0,100,0
250     100,0,0
@HuidaeCho

This comment has been minimized.

Copy link

HuidaeCho commented May 17, 2014

Combined with ticket #4918, the following patch fixes this issue:

--- mapdrawgdal.c.orig  2014-05-16 05:00:53.811364470 -0400
+++ mapdrawgdal.c   2014-05-16 06:25:48.419077829 -0400
@@ -49,13 +49,18 @@
 #include "gdal_alg.h"

 static int
+IsNoData( double dfValue, double dfNoDataValue, int bIsNoDataNan );
+
+static int
 LoadGDALImages( GDALDatasetH hDS, int band_numbers[4], int band_count,
                 layerObj *layer,
                 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
                 GByte *pabyBuffer,
                 int dst_xsize, int dst_ysize,
                 int *pbHaveRGBNoData,
-                int *pnNoData1, int *pnNoData2, int *pnNoData3 );
+                int *pnNoData1, int *pnNoData2, int *pnNoData3,
+                int *scaled, double *scaleMin, double *scaleRatio,
+                unsigned char **isNoDataBuffer );
 static int
 msDrawRasterLayerGDAL_RawMode(
   mapObj *map, layerObj *layer, imageObj *image, GDALDatasetH hDS,
@@ -131,6 +136,11 @@
   int bHaveRGBNoData = FALSE;
   int nNoData1=-1,nNoData2=-1,nNoData3=-1;
   rasterBufferObj *mask_rb = NULL;
+  int isDefaultGreyscale;
+  int scaled;
+  double scaleMin;
+  double scaleRatio;
+  unsigned char *isNoDataBuffer;
 #ifdef USE_GD
   int   anColorCube[256];
   int cmt=0;
@@ -492,6 +502,8 @@
    * Get colormap for this image.  If there isn't one, and we have only
    * one band create a greyscale colormap.
    */
+  isDefaultGreyscale = FALSE;
+
   if( hBand2 != NULL )
     hColorMap = NULL;
   else {
@@ -500,6 +512,7 @@
       hColorMap = GDALCloneColorTable( hColorMap );
     else if( hBand2 == NULL ) {
       hColorMap = GDALCreateColorTable( GPI_RGB );
+      isDefaultGreyscale = TRUE;

       for( i = 0; i < 256; i++ ) {
         colorObj pixel;
@@ -561,13 +574,49 @@
   }

   /*
+   * Allocate imagery buffers.
+   */
+  pabyRaw1 = (unsigned char *) malloc(dst_xsize * dst_ysize * band_count);
+  if( pabyRaw1 == NULL ) {
+    msSetError(MS_MEMERR, "Allocating work image of size %dx%dx%d failed.",
+               "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize, band_count );
+    return -1;
+  }
+
+  if( hBand2 != NULL && hBand3 != NULL ) {
+    pabyRaw2 = pabyRaw1 + dst_xsize * dst_ysize * 1;
+    pabyRaw3 = pabyRaw1 + dst_xsize * dst_ysize * 2;
+  }
+
+  if( hBandAlpha != NULL ) {
+    if( hBand2 != NULL )
+      pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 3;
+    else
+      pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 1;
+  }
+
+  /*
+   * Load image data into buffers with scaling, etc.
+   */
+  if( LoadGDALImages( hDS, band_numbers, band_count, layer,
+                      src_xoff, src_yoff, src_xsize, src_ysize,
+                      pabyRaw1, dst_xsize, dst_ysize,
+                      &bHaveRGBNoData,
+                      &nNoData1, &nNoData2, &nNoData3,
+             &scaled, &scaleMin, &scaleRatio,
+             &isNoDataBuffer ) == -1 ) {
+    free( pabyRaw1 );
+    return -1;
+  }
+
+  /*
    * Setup the mapping between source eight bit pixel values, and the
    * output images color table.  There are two general cases, where the
    * class colors are provided by the MAP file, or where we use the native
    * color table.
    */
   if( classified ) {
-    int c, color_count;
+    int c, color_count, scaleColors;

 #ifndef NDEBUG
     cmap_set = TRUE;
@@ -580,7 +629,12 @@
       return -1;
     }

-    color_count = MIN(256,GDALGetColorEntryCount(hColorMap));
+    color_count = GDALGetColorEntryCount(hColorMap);
+    scaleColors = scaled && !isDefaultGreyscale;
+
+    if( !scaleColors && color_count > 256 )
+      color_count = 256;
+
     for(i=0; i < color_count; i++) {
       colorObj pixel;
       int colormap_index;
@@ -588,17 +642,19 @@

       GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );

+      j = scaleColors ? MAX(0, MIN(255, (int) ((i-scaleMin)*scaleRatio))) : i;
+
       pixel.red = sEntry.c1;
       pixel.green = sEntry.c2;
       pixel.blue = sEntry.c3;
-      colormap_index = i;
+      colormap_index = j;

       if(!MS_COMPARE_COLORS(pixel, layer->offsite)) {
         c = msGetClass(layer, &pixel, colormap_index);

         if(c == -1) { /* doesn't belong to any class, so handle like offsite*/
           if( rb->type == MS_BUFFER_GD )
-            cmap[i] = -1;
+            cmap[j] = -1;
         } else {
           int s;

@@ -615,12 +671,12 @@
           if( rb->type == MS_BUFFER_GD ) {
             RESOLVE_PEN_GD(rb->data.gd_img, layer->class[c]->styles[0]->color);
             if( MS_TRANSPARENT_COLOR(layer->class[c]->styles[0]->color) )
-              cmap[i] = -1;
+              cmap[j] = -1;
             else if( MS_VALID_COLOR(layer->class[c]->styles[0]->color)) {
               /* use class color */
-              cmap[i] = layer->class[c]->styles[0]->color.pen;
+              cmap[j] = layer->class[c]->styles[0]->color.pen;
             } else /* Use raster color */
-              cmap[i] = msAddColorGD(map, rb->data.gd_img, cmt,
+              cmap[j] = msAddColorGD(map, rb->data.gd_img, cmt,
                                      pixel.red, pixel.green, pixel.blue);
           } else if( rb->type == MS_BUFFER_BYTE_RGBA )
 #endif
@@ -629,74 +685,86 @@
               /* leave it transparent */;

             else if( MS_VALID_COLOR(layer->class[c]->styles[0]->color)) {
-              rb_cmap[0][i] = layer->class[c]->styles[0]->color.red;
-              rb_cmap[1][i] = layer->class[c]->styles[0]->color.green;
-              rb_cmap[2][i] = layer->class[c]->styles[0]->color.blue;
-              rb_cmap[3][i] = (255*layer->class[c]->styles[0]->opacity / 100);
+              rb_cmap[0][j] = layer->class[c]->styles[0]->color.red;
+              rb_cmap[1][j] = layer->class[c]->styles[0]->color.green;
+              rb_cmap[2][j] = layer->class[c]->styles[0]->color.blue;
+              rb_cmap[3][j] = (255*layer->class[c]->styles[0]->opacity / 100);
             }

             else { /* Use raster color */
-              rb_cmap[0][i] = pixel.red;
-              rb_cmap[1][i] = pixel.green;
-              rb_cmap[2][i] = pixel.blue;
-              rb_cmap[3][i] = 255;
+              rb_cmap[0][j] = pixel.red;
+              rb_cmap[1][j] = pixel.green;
+              rb_cmap[2][j] = pixel.blue;
+              rb_cmap[3][j] = 255;
             }
           }
         }
 #ifdef USE_GD
       } else {
         if( rb->type == MS_BUFFER_GD )
-          cmap[i] = -1;
+          cmap[j] = -1;
 #endif
       }
     }
 #ifdef USE_GD
   } else if( hColorMap != NULL && rb->type == MS_BUFFER_GD ) {
-    int color_count;
+    int color_count, scaleColors;
 #ifndef NDEBUG
     cmap_set = TRUE;
 #endif

-    color_count = MIN(256,GDALGetColorEntryCount(hColorMap));
+    color_count = GDALGetColorEntryCount(hColorMap);
+    scaleColors = scaled && !isDefaultGreyscale;
+
+    if( !scaleColors && color_count > 256 )
+      color_count = 256;

     for(i=0; i < color_count; i++) {
       GDALColorEntry sEntry;

       GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );

+      j = scaleColors ? MAX(0, MIN(255, (int) ((i-scaleMin)*scaleRatio))) : i;
+
       if( sEntry.c4 != 0
           && (!MS_VALID_COLOR( layer->offsite )
               || layer->offsite.red != sEntry.c1
               || layer->offsite.green != sEntry.c2
               || layer->offsite.blue != sEntry.c3 ) )
-        cmap[i] = msAddColorGD(map, rb->data.gd_img, cmt,
+        cmap[j] = msAddColorGD(map, rb->data.gd_img, cmt,
                                sEntry.c1, sEntry.c2, sEntry.c3);
       else
-        cmap[i] = -1;
+        cmap[j] = -1;
     }
 #endif
   } else if( hBand2 == NULL && hColorMap != NULL && rb->type == MS_BUFFER_BYTE_RGBA ) {
-    int color_count;
+    int color_count, scaleColors;
 #ifndef NDEBUG
     cmap_set = TRUE;
 #endif

-    color_count = MIN(256,GDALGetColorEntryCount(hColorMap));
+    color_count = GDALGetColorEntryCount(hColorMap);
+    scaleColors = scaled && !isDefaultGreyscale;
+
+    if( !scaleColors && color_count > 256 )
+      color_count = 256;

     for(i=0; i < color_count; i++) {
       GDALColorEntry sEntry;

       GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );

+      j = scaleColors ? MAX(0, MIN(255, (int) ((i-scaleMin)*scaleRatio))) : i;
+
       if( sEntry.c4 != 0
           && (!MS_VALID_COLOR( layer->offsite )
               || layer->offsite.red != sEntry.c1
               || layer->offsite.green != sEntry.c2
               || layer->offsite.blue != sEntry.c3 ) ) {
-        rb_cmap[0][i] = sEntry.c1;
-        rb_cmap[1][i] = sEntry.c2;
-        rb_cmap[2][i] = sEntry.c3;
-        rb_cmap[3][i] = sEntry.c4;
+        rb_cmap[0][j] = sEntry.c1;
+        rb_cmap[1][j] = sEntry.c2;
+        rb_cmap[2][j] = sEntry.c3;
+        rb_cmap[3][j] = sEntry.c4;
       }
     }
   }
@@ -706,40 +774,6 @@
   }
 #endif

-  /*
-   * Allocate imagery buffers.
-   */
-  pabyRaw1 = (unsigned char *) malloc(dst_xsize * dst_ysize * band_count);
-  if( pabyRaw1 == NULL ) {
-    msSetError(MS_MEMERR, "Allocating work image of size %dx%dx%d failed.",
-               "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize, band_count );
-    return -1;
-  }
-
-  if( hBand2 != NULL && hBand3 != NULL ) {
-    pabyRaw2 = pabyRaw1 + dst_xsize * dst_ysize * 1;
-    pabyRaw3 = pabyRaw1 + dst_xsize * dst_ysize * 2;
-  }
-
-  if( hBandAlpha != NULL ) {
-    if( hBand2 != NULL )
-      pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 3;
-    else
-      pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 1;
-  }
-
-  /*
-   * Load image data into buffers with scaling, etc.
-   */
-  if( LoadGDALImages( hDS, band_numbers, band_count, layer,
-                      src_xoff, src_yoff, src_xsize, src_ysize,
-                      pabyRaw1, dst_xsize, dst_ysize,
-                      &bHaveRGBNoData,
-                      &nNoData1, &nNoData2, &nNoData3 ) == -1 ) {
-    free( pabyRaw1 );
-    return -1;
-  }
-
   if( bHaveRGBNoData && layer->debug )
     msDebug( "msDrawGDAL(): using RGB nodata values from GDAL dataset.\n" );

@@ -868,12 +902,10 @@

       k = 0;
       for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
-        for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
+        for( j = dst_xoff; j < dst_xoff + dst_xsize; j++, k++ ) {
           int src_pixel, src_alpha, cmap_alpha, merged_alpha;
-          if(SKIP_MASK(j,i)) {
-            k++;
+          if(SKIP_MASK(j,i) || (isNoDataBuffer && isNoDataBuffer[k]))
             continue;
-          }

           src_pixel = pabyRaw1[k];
           src_alpha = pabyRawAlpha[k];
@@ -896,7 +928,6 @@
                           rb_cmap[2][src_pixel],
                           merged_alpha );
           }
-          k++;
         }
       }
     }
@@ -909,11 +940,10 @@

       k = 0;
       for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
-        for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
-          int src_pixel = pabyRaw1[k++];
-          if(SKIP_MASK(j,i)) {
+        for( j = dst_xoff; j < dst_xoff + dst_xsize; j++, k++ ) {
+          int src_pixel = pabyRaw1[k];
+          if(SKIP_MASK(j,i) || (isNoDataBuffer && isNoDataBuffer[k]))
             continue;
-          }

           if( rb_cmap[3][src_pixel] > 253 ) {
             RB_SET_PIXEL( rb, j, i,
@@ -1066,6 +1096,9 @@
   msFree( mask_rb );
   free( pabyRaw1 );

+  if( isNoDataBuffer )
+    free( isNoDataBuffer );
+
   if( hColorMap != NULL )
     GDALDestroyColorTable( hColorMap );

@@ -1329,13 +1362,17 @@
                 GByte *pabyWholeBuffer,
                 int dst_xsize, int dst_ysize,
                 int *pbHaveRGBNoData,
-                int *pnNoData1, int *pnNoData2, int *pnNoData3 )
-
+                int *pnNoData1, int *pnNoData2, int *pnNoData3,
+                int *scaled, double *scaleMin, double *scaleRatio,
+                unsigned char **isNoDataBuffer)
 {
   int    iColorIndex, result_code=0;
   CPLErr eErr;
   float *pafWholeRawData;

+  *scaled = FALSE;
+  *isNoDataBuffer = NULL;
+
   /* -------------------------------------------------------------------- */
   /*      If we have no alpha band, but we do have three input            */
   /*      bands, then check for nodata values.  If we only have one       */
@@ -1446,9 +1483,11 @@
   for( iColorIndex = 0; iColorIndex < band_count; iColorIndex++ ) {
     unsigned char *pabyBuffer;
     double dfScaleMin=0.0, dfScaleMax=255.0, dfScaleRatio, dfNoDataValue;
+    float fNoDataValue;
     const char *pszScaleInfo;
     float *pafRawData;
     int   nPixelCount = dst_xsize * dst_ysize, i, bGotNoData = FALSE;
+    int bIsNoDataNan;
     GDALRasterBandH hBand =GDALGetRasterBand(hDS,band_numbers[iColorIndex]);
     pszScaleInfo = CSLFetchNameValue( layer->processing, "SCALE" );
     if( pszScaleInfo == NULL ) {
@@ -1489,16 +1528,17 @@
     pafRawData = pafWholeRawData + iColorIndex * dst_xsize * dst_ysize;

     dfNoDataValue = msGetGDALNoDataValue( layer, hBand, &bGotNoData );
+    bIsNoDataNan = isnan(dfNoDataValue);
+
+    /* we force assignment to a float rather than letting pafRawData[i]
+       get promoted to double later to avoid float precision issues. */
+    fNoDataValue = (float) dfNoDataValue;

     if( dfScaleMin == dfScaleMax ) {
       int bMinMaxSet = 0;

-      /* we force assignment to a float rather than letting pafRawData[i]
-         get promoted to double later to avoid float precision issues. */
-      float fNoDataValue = (float) dfNoDataValue;
-
       for( i = 0; i < nPixelCount; i++ ) {
-        if( bGotNoData && pafRawData[i] == fNoDataValue )
+        if( bGotNoData && IsNoData(pafRawData[i], fNoDataValue, bIsNoDataNan) )
           continue;

         if( !bMinMaxSet ) {
@@ -1524,8 +1564,19 @@
     dfScaleRatio = 256.0 / (dfScaleMax - dfScaleMin);
     pabyBuffer = pabyWholeBuffer + iColorIndex * nPixelCount;

+    if( iColorIndex == 0 && bGotNoData )
+      *isNoDataBuffer = (unsigned char *)calloc(nPixelCount, 1);
+
     for( i = 0; i < nPixelCount; i++ ) {
-      float fScaledValue = (float) ((pafRawData[i]-dfScaleMin)*dfScaleRatio);
+      float fScaledValue;
+
+      if( bGotNoData && IsNoData(pafRawData[i], fNoDataValue, bIsNoDataNan) ) {
+   if( iColorIndex == 0 )
+     (*isNoDataBuffer)[i] = 1;
+        continue;
+      }
+
+      fScaledValue = (float) ((pafRawData[i]-dfScaleMin)*dfScaleRatio);

       if( fScaledValue < 0.0 )
         pabyBuffer[i] = 0;
@@ -1535,6 +1586,12 @@
         pabyBuffer[i] = (int) fScaledValue;
     }

+    if( iColorIndex == 0 ) {
+      *scaled = TRUE;
+      *scaleMin = dfScaleMin;
+      *scaleRatio = dfScaleRatio;
+    }
+
     /* -------------------------------------------------------------------- */
     /*      Report a warning if NODATA keyword was applied.  We are         */
     /*      unable to utilize it since we can't return any pixels marked    */
@@ -2088,6 +2145,7 @@
   const char *pszScaleInfo;
   const char *pszBuckets;
   int  *cmap, c, j, k, bGotNoData = FALSE, bGotFirstValue;
+  int bIsNoDataNan;
   unsigned char *rb_cmap[4];
   CPLErr eErr;
   rasterBufferObj *mask_rb = NULL;
@@ -2126,6 +2184,7 @@
   }

   fNoDataValue = (float) msGetGDALNoDataValue( layer, hBand, &bGotNoData );
+  bIsNoDataNan = isnan(fNoDataValue);

   /* ==================================================================== */
   /*      Determine scaling.                                              */
@@ -2138,7 +2197,7 @@
   bGotFirstValue = FALSE;

   for( i = 1; i < nPixelCount; i++ ) {
-    if( bGotNoData && pafRawData[i] == fNoDataValue )
+    if( bGotNoData && IsNoData(pafRawData[i], fNoDataValue, bIsNoDataNan) )
       continue;

     if( !bGotFirstValue ) {
@@ -2301,7 +2360,7 @@
       /*
        * Skip nodata pixels ... no processing.
        */
-      if( bGotNoData && fRawValue == fNoDataValue ) {
+      if( bGotNoData && IsNoData(fRawValue, fNoDataValue, bIsNoDataNan) ) {
         continue;
       }

@@ -2355,6 +2414,16 @@
 }

 /************************************************************************/
+/*                          IsNoData()                                  */
+/************************************************************************/
+static int
+IsNoData( double dfValue, double dfNoDataValue, int bIsNoDataNan )
+{
+  return (bIsNoDataNan && isnan(dfValue)) ||
+    (!bIsNoDataNan && dfValue == dfNoDataValue);
+}
+
+/************************************************************************/
 /*                          msGetGDALNoDataValue()                      */
 /************************************************************************/
 double
@sdlime

This comment has been minimized.

Copy link
Member

sdlime commented May 20, 2014

@warmerdam, any chance you could review quickly? --Steve

@mapserver-bot

This comment has been minimized.

Copy link

mapserver-bot commented Feb 23, 2016

This is an automated comment

This issue has been closed due to lack of activity. This doesn't mean the issue is invalid, it simply got no attention within the last year. Please reopen with missing/relevant information if still valid.

Typically, issues fall in this state for one of the following reasons:

  • Hard, impossible or not enough information to reproduce
  • Missing test case
  • Lack of a champion with interest and/or funding to address the issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment