root/tags/rel-5-0-2/mapserver/mapdrawgdal.c

/******************************************************************************
 *
 * Project:  MapServer
 * Purpose:  Code for drawing GDAL raster layers.  Called from 
 *           msDrawRasterLayerLow() in mapraster.c.  
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in 
 * all copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#include <assert.h>
#include "mapserver.h"
#include "mapresample.h"
#include "mapthread.h"

MS_CVSID("$Id$")

extern int InvGeoTransform( double *gt_in, double *gt_out );

#define MAXCOLORS 256
#define GEO_TRANS(tr,x,y)  ((tr)[0]+(tr)[1]*(x)+(tr)[2]*(y))

#if defined(USE_GDAL)

#include "gdal.h"
#include "cpl_string.h"

#if GDAL_VERSION_NUM > 1174
#  define ENABLE_DITHER
#endif

#ifdef ENABLE_DITHER
#include "gdal_alg.h"
#endif

static int
LoadGDALImage( GDALRasterBandH hBand, int iColorIndex, 
               layerObj *layer, 
               int src_xoff, int src_yoff, int src_xsize, int src_ysize, 
               GByte *pabyBuffer, int dst_xsize, int dst_ysize );
static int 
msDrawRasterLayerGDAL_RawMode(
    mapObj *map, layerObj *layer, imageObj *image, GDALDatasetH hDS, 
    int src_xoff, int src_yoff, int src_xsize, int src_ysize,
    int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize );

static int 
msDrawRasterLayerGDAL_16BitClassification(
    mapObj *map, layerObj *layer, imageObj *image, 
    GDALDatasetH hDS, GDALRasterBandH hBand,
    int src_xoff, int src_yoff, int src_xsize, int src_ysize,
    int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize );


#ifdef ENABLE_DITHER
static void Dither24to8( GByte *pabyRed, GByte *pabyGreen, GByte *pabyBlue,
                         GByte *pabyDithered, int xsize, int ysize, 
                         int bTransparent, colorObj transparentColor,
                         gdImagePtr gdImg );
#endif

/*
 * Stuff for allocating color cubes, and mapping between RGB values and
 * the corresponding color cube index.
 */

static int allocColorCube(mapObj *map, gdImagePtr img, int *panColorCube);

#define RED_LEVELS 5
#define RED_DIV 52
/* #define GREEN_LEVELS 7 */
#define GREEN_LEVELS 5
/* #define GREEN_DIV 37 */
#define GREEN_DIV 52
#define BLUE_LEVELS 5
#define BLUE_DIV 52

#define RGB_LEVEL_INDEX(r,g,b) ((r)*GREEN_LEVELS*BLUE_LEVELS + (g)*BLUE_LEVELS+(b))
#define RGB_INDEX(r,g,b) RGB_LEVEL_INDEX(((r)/RED_DIV),((g)/GREEN_DIV),((b)/BLUE_DIV))

/************************************************************************/
/*                       msDrawRasterLayerGDAL()                        */
/************************************************************************/

int msDrawRasterLayerGDAL(mapObj *map, layerObj *layer, imageObj *image, 
                          void *hDSVoid )

{
  int i,j, k; /* loop counters */
  int cmap[MAXCOLORS], cmap_set = FALSE;
  double adfGeoTransform[6], adfInvGeoTransform[6];
  int   dst_xoff, dst_yoff, dst_xsize, dst_ysize;
  int   src_xoff, src_yoff, src_xsize, src_ysize;
  int   anColorCube[256];
  double llx, lly, urx, ury;
  rectObj copyRect, mapRect;
  unsigned char *pabyRaw1=NULL, *pabyRaw2=NULL, *pabyRaw3=NULL, 
                *pabyRawAlpha = NULL;
  int truecolor = FALSE, classified = FALSE;
  int red_band=0, green_band=0, blue_band=0, alpha_band=0, cmt=0;
  gdImagePtr gdImg = NULL;
  GDALDatasetH hDS = hDSVoid;
  GDALColorTableH hColorMap;
  GDALRasterBandH hBand1, hBand2, hBand3, hBandAlpha;

  memset( cmap, 0xff, MAXCOLORS * sizeof(int) );

/* -------------------------------------------------------------------- */
/*      Test the image format instead of the map format.                */
/*      Normally the map format and the image format should be the      */
/*      same but In somes cases like swf and pdf support, a temporary   */
/*      GD image object is created and used to render raster layers     */
/*      and then dumped into the SWF or the PDF file.                   */
/* -------------------------------------------------------------------- */
  
  /* if( MS_RENDERER_GD(map->outputformat) ) */
  if( MS_RENDERER_GD(image->format) )
  {
      gdImg = image->img.gd;

      truecolor = gdImageTrueColor( gdImg );
      if( CSLFetchNameValue( layer->processing, 
                             "COLOR_MATCH_THRESHOLD" ) != NULL )
      {
          cmt = MAX(0,atoi(CSLFetchNameValue( layer->processing, 
                                              "COLOR_MATCH_THRESHOLD" )));
      }
  }
#ifdef USE_AGG
  else if( MS_RENDERER_AGG(image->format) )
  {
      gdImg = image->img.gd;

      truecolor = gdImageTrueColor( gdImg );
      if( CSLFetchNameValue( layer->processing, 
                             "COLOR_MATCH_THRESHOLD" ) != NULL )
      {
          cmt = MAX(0,atoi(CSLFetchNameValue( layer->processing, 
                                              "COLOR_MATCH_THRESHOLD" )));
      }
  }
#endif
  src_xsize = GDALGetRasterXSize( hDS );
  src_ysize = GDALGetRasterYSize( hDS );

  /*
   * If the RAW_WINDOW attribute is set, use that to establish what to
   * load.  This is normally just set by the mapresample.c module to avoid
   * problems with rotated maps.
   */

  if( CSLFetchNameValue( layer->processing, "RAW_WINDOW" ) != NULL )
  {
      char **papszTokens = 
          CSLTokenizeString( 
              CSLFetchNameValue( layer->processing, "RAW_WINDOW" ) );
      
      if( layer->debug )
          msDebug( "msDrawGDAL(%s): using RAW_WINDOW=%s\n",
                   layer->name,
                   CSLFetchNameValue( layer->processing, "RAW_WINDOW" ) );

      if( CSLCount(papszTokens) != 4 )
      {
          CSLDestroy( papszTokens );
          msSetError( MS_IMGERR, "RAW_WINDOW PROCESSING directive corrupt.",
                      "msDrawGDAL()" );
          return -1;
      }

      src_xoff = atoi(papszTokens[0]);
      src_yoff = atoi(papszTokens[1]);
      src_xsize = atoi(papszTokens[2]);
      src_ysize = atoi(papszTokens[3]);

      dst_xoff = 0;
      dst_yoff = 0;
      dst_xsize = image->width;
      dst_ysize = image->height;

      CSLDestroy( papszTokens );
  }

  /*
   * Compute the georeferenced window of overlap, and do nothing if there
   * is no overlap between the map extents, and the file we are operating on.
   */
  else if( layer->transform )
  {
      int dst_lrx, dst_lry;

      msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
      InvGeoTransform( adfGeoTransform, adfInvGeoTransform );
      
      mapRect = map->extent;
      
      mapRect.minx -= map->cellsize*0.5;
      mapRect.maxx += map->cellsize*0.5;
      mapRect.miny -= map->cellsize*0.5;
      mapRect.maxy += map->cellsize*0.5;
      
      copyRect = mapRect;
      
      if( copyRect.minx < GEO_TRANS(adfGeoTransform,0,src_ysize) )
          copyRect.minx = GEO_TRANS(adfGeoTransform,0,src_ysize);
      if( copyRect.maxx > GEO_TRANS(adfGeoTransform,src_xsize,0) )
          copyRect.maxx = GEO_TRANS(adfGeoTransform,src_xsize,0);

      if( copyRect.miny < GEO_TRANS(adfGeoTransform+3,0,src_ysize) )
          copyRect.miny = GEO_TRANS(adfGeoTransform+3,0,src_ysize);
      if( copyRect.maxy > GEO_TRANS(adfGeoTransform+3,src_xsize,0) )
          copyRect.maxy = GEO_TRANS(adfGeoTransform+3,src_xsize,0);
      
      if( copyRect.minx >= copyRect.maxx || copyRect.miny >= copyRect.maxy )
          return 0;

      /*
       * Copy the source and destination raster coordinates.
       */
      llx = GEO_TRANS(adfInvGeoTransform+0,copyRect.minx,copyRect.miny);
      lly = GEO_TRANS(adfInvGeoTransform+3,copyRect.minx,copyRect.miny);
      urx = GEO_TRANS(adfInvGeoTransform+0,copyRect.maxx,copyRect.maxy);
      ury = GEO_TRANS(adfInvGeoTransform+3,copyRect.maxx,copyRect.maxy);
      
      src_xoff = MAX(0,(int) llx);
      src_yoff = MAX(0,(int) ury);
      src_xsize = MIN(MAX(0,(int) (urx - llx + 0.5)),
                      GDALGetRasterXSize(hDS) - src_xoff);
      src_ysize = MIN(MAX(0,(int) (lly - ury + 0.5)),
                      GDALGetRasterYSize(hDS) - src_yoff);

      if( src_xsize == 0 || src_ysize == 0 )
      {
          if( layer->debug )
              msDebug( "msDrawGDAL(): no apparent overlap between map view and this window(1).\n" );
          return 0;
      }

      dst_xoff = (int) ((copyRect.minx - mapRect.minx) / map->cellsize);
      dst_yoff = (int) ((mapRect.maxy - copyRect.maxy) / map->cellsize);

      dst_lrx = (int) ((copyRect.maxx - mapRect.minx) / map->cellsize + 0.5);
      dst_lry = (int) ((mapRect.maxy - copyRect.miny) / map->cellsize + 0.5);
      dst_lrx = MAX(0,MIN(image->width,dst_lrx));
      dst_lry = MAX(0,MIN(image->height,dst_lry));
      
      dst_xsize = MAX(0,MIN(image->width,dst_lrx - dst_xoff));
      dst_ysize = MAX(0,MIN(image->height,dst_lry - dst_yoff));

      if( dst_xsize == 0 || dst_ysize == 0 )
      {
          if( layer->debug )
              msDebug( "msDrawGDAL(): no apparent overlap between map view and this window(2).\n" );
          return 0;
      }

#ifndef notdef
      if( layer->debug )
          msDebug( "msDrawGDAL(): src=%d,%d,%d,%d, dst=%d,%d,%d,%d\n", 
                   src_xoff, src_yoff, src_xsize, src_ysize, 
                   dst_xoff, dst_yoff, dst_xsize, dst_ysize );
#endif
  }

  /*
   * If layer transforms are turned off, just map 1:1.
   */
  else
  {
      dst_xoff = src_xoff = 0;
      dst_yoff = src_yoff = 0;
      dst_xsize = src_xsize = MIN(image->width,src_xsize);
      dst_ysize = src_ysize = MIN(image->height,src_ysize);
  }

  /*
   * In RAWDATA mode we don't fool with colors.  Do the raw processing, 
   * and return from the function early.
   */
  if( !gdImg )
  {
      assert( MS_RENDERER_RAWDATA( image->format ) );

      return msDrawRasterLayerGDAL_RawMode( 
          map, layer, image, hDS, 
          src_xoff, src_yoff, src_xsize, src_ysize, 
          dst_xoff, dst_yoff, dst_xsize, dst_ysize );
  }
  
  /*
   * Is this image classified?  We consider it classified if there are
   * classes with an expression string *or* a color range.  We don't want
   * to treat the raster as classified if there is just a bogus class here
   * to force inclusion in the legend.
   */
  for( i = 0; i < layer->numclasses; i++ )
  {
      int s;

      /* change colour based on colour range? */
      for(s=0; s<layer->class[i]->numstyles; s++)
      {
          if( MS_VALID_COLOR(layer->class[i]->styles[s]->mincolor)
              && MS_VALID_COLOR(layer->class[i]->styles[s]->maxcolor) )
          {
              classified = TRUE;
              break;
          }
      }
      
      if( layer->class[i]->expression.string != NULL )
      {
          classified = TRUE;
          break;
      }
  }

  /*
   * Set up the band selection.  We look for a BANDS directive in the 
   * the PROCESSING options.  If not found we default to red=1 or
   * red=1,green=2,blue=3 or red=1,green=2,blue=3,alpha=4. 
   */

  if( CSLFetchNameValue( layer->processing, "BANDS" ) == NULL )
  {
      red_band = 1;

      if( GDALGetRasterCount( hDS ) >= 4 
          && GDALGetRasterColorInterpretation( 
              GDALGetRasterBand( hDS, 4 ) ) == GCI_AlphaBand )
          alpha_band = 4;
      
      if( GDALGetRasterCount( hDS ) >= 3 )
      {
          green_band = 2;
          blue_band = 3;
      }

      if( GDALGetRasterCount( hDS ) == 2 
          && GDALGetRasterColorInterpretation( 
              GDALGetRasterBand( hDS, 2 ) ) == GCI_AlphaBand )
          alpha_band = 2;
      
      hBand1 = GDALGetRasterBand( hDS, red_band );
      if( classified 
          || GDALGetRasterColorTable( hBand1 ) != NULL )
      {
          alpha_band = 0;
          green_band = 0;
          blue_band = 0;
      }
  }
  else
  {
      int band_count, *band_list;

      band_list = msGetGDALBandList( layer, hDS, 4, &band_count );
      if( band_list == NULL )
          return -1;
      
      if( band_count > 0 )
          red_band = band_list[0];
      else
          red_band = 0;
      if( band_count > 2 )
      {
          green_band = band_list[1];
          blue_band = band_list[2];
      }
      else      
      {
          green_band = 0;
          blue_band = 0;
      }

      if( band_count > 3 )
          alpha_band = band_list[3];
      else
          alpha_band = 0;

      free( band_list );
  }

  if( layer->debug > 1 || (layer->debug > 0 && green_band != 0) )
  {
      msDebug( "msDrawGDAL(): red,green,blue,alpha bands = %d,%d,%d,%d\n", 
               red_band, green_band, blue_band, alpha_band );
  }

  /*
   * Get band handles for PC256, RGB or RGBA cases.
   */
  hBand1 = GDALGetRasterBand( hDS, red_band );
  if( hBand1 == NULL )
      return -1;

  hBand2 = hBand3 = hBandAlpha = NULL;

  if( green_band != 0 )
  {
      hBand1 = GDALGetRasterBand( hDS, red_band );
      hBand2 = GDALGetRasterBand( hDS, green_band );
      hBand3 = GDALGetRasterBand( hDS, blue_band );
      if( hBand1 == NULL || hBand2 == NULL || hBand3 == NULL )
          return -1;
  }

  if( alpha_band != 0 )
      hBandAlpha = GDALGetRasterBand( hDS, alpha_band );

  /*
   * Wipe pen indicators for all our layer class colors if they exist.  
   * Sometimes temporary gdImg'es are used in which case previously allocated
   * pens won't generally apply.  See Bug 504.
   */
  if( gdImg && !truecolor )
  {
      int iClass;
      int iStyle;
      for( iClass = 0; iClass < layer->numclasses; iClass++ )
      {
          for (iStyle=0; iStyle<layer->class[iClass]->numstyles; iStyle++)
            layer->class[iClass]->styles[iStyle]->color.pen = MS_PEN_UNSET;
      }
  }

  /*
   * The logic for a classification rendering of non-8bit raster bands
   * is sufficiently different than the normal mechanism of loading
   * into an 8bit buffer, that we isolate it into it's own subfunction.
   */
  if( classified && gdImg 
     && hBand1 != NULL && GDALGetRasterDataType( hBand1 ) != GDT_Byte ) 
  {
      return msDrawRasterLayerGDAL_16BitClassification( 
          map, layer, image, hDS, hBand1,
          src_xoff, src_yoff, src_xsize, src_ysize, 
          dst_xoff, dst_yoff, dst_xsize, dst_ysize );
  }

  /*
   * Get colormap for this image.  If there isn't one, and we have only
   * one band create a greyscale colormap. 
   */
  if( hBand2 != NULL )
      hColorMap = NULL;
  else
  {
      hColorMap = GDALGetRasterColorTable( hBand1 );
      if( hColorMap != NULL )
          hColorMap = GDALCloneColorTable( hColorMap );
      else if( hBand2 == NULL )
      {
          hColorMap = GDALCreateColorTable( GPI_RGB );
          
          for( i = 0; i < 256; i++ )
          {
              colorObj pixel;
              GDALColorEntry sEntry;

              pixel.red = i;
              pixel.green = i;
              pixel.blue = i;
              pixel.pen = i;
              
              if(MS_COMPARE_COLORS(pixel, layer->offsite))
              {
                  sEntry.c1 = 0;
                  sEntry.c2 = 0;
                  sEntry.c3 = 0;
                  sEntry.c4 = 0; /* alpha set to zero */
              }
              else if( truecolor )
              {
                  sEntry.c1 = i;
                  sEntry.c2 = i;
                  sEntry.c3 = i;
                  sEntry.c4 = 255;
              }
              else
              {
                  /*
                  ** This special calculation is intended to use only 128
                  ** unique colors for greyscale in non-truecolor mode.
                  */

                  sEntry.c1 = i - i%2;
                  sEntry.c2 = i - i%2;
                  sEntry.c3 = i - i%2;
                  sEntry.c4 = 255;
              }
                  
              GDALSetColorEntry( hColorMap, i, &sEntry );
          }
      }

      /* 
      ** If we have a known NODATA value, mark it now as transparent. 
      */
      {
          int    bGotNoData;
          double dfNoDataValue = msGetGDALNoDataValue( layer, hBand1, 
                                                       &bGotNoData);

          if( bGotNoData && dfNoDataValue >= 0 
              && dfNoDataValue < GDALGetColorEntryCount( hColorMap ) )
          {
              GDALColorEntry sEntry;
              
              memcpy( &sEntry, 
                      GDALGetColorEntry( hColorMap, (int) dfNoDataValue ),
                      sizeof(GDALColorEntry) );

              sEntry.c4 = 0;
              GDALSetColorEntry( hColorMap, (int) dfNoDataValue, &sEntry );
          }
      }
  }

  /*
   * Setup the mapping between source eightbit 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 && gdImg ) {
    int c, color_count;

    cmap_set = TRUE;

    if( hColorMap == NULL )
    {
        msSetError(MS_IOERR, 
                   "Attempt to classify 24bit image, this is unsupported.",
                   "drawGDAL()");
        return -1;
    }

    color_count = MIN(256,GDALGetColorEntryCount(hColorMap));
    for(i=0; i < color_count; i++) {
        colorObj pixel;
        GDALColorEntry sEntry;

        GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );
            
        pixel.red = sEntry.c1;
        pixel.green = sEntry.c2;
        pixel.blue = sEntry.c3;
        pixel.pen = i;
        
        if(!MS_COMPARE_COLORS(pixel, layer->offsite))
        {
            c = msGetClass(layer, &pixel);
            
            if(c == -1)/* doesn't belong to any class, so handle like offsite*/
                cmap[i] = -1;
            else
            {
                int s;
                
                /* change colour based on colour range?  Currently we 
                   only address the greyscale case properly. */

                for(s=0; s<layer->class[c]->numstyles; s++)
                {
                    if( MS_VALID_COLOR(layer->class[c]->styles[s]->mincolor)
                        && MS_VALID_COLOR(layer->class[c]->styles[s]->maxcolor) )
                        msValueToRange(layer->class[c]->styles[s],
                                       sEntry.c1 );
                }

                RESOLVE_PEN_GD(gdImg, layer->class[c]->styles[0]->color);
                if( MS_TRANSPARENT_COLOR(layer->class[c]->styles[0]->color) )
                    cmap[i] = -1;
                else if( MS_VALID_COLOR(layer->class[c]->styles[0]->color))
                {
                    /* use class color */
                    cmap[i] = layer->class[c]->styles[0]->color.pen;
                }
                else /* Use raster color */
                    cmap[i] = msAddColorGD(map, gdImg, cmt,
                                           pixel.red, pixel.green, pixel.blue);
            }
        } else
            cmap[i] = -1;
    }
  } else if( hColorMap != NULL && !truecolor && gdImg ) {
    int color_count;
    cmap_set = TRUE;

    color_count = MIN(256,GDALGetColorEntryCount(hColorMap));

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

        GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );

        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, gdImg, cmt, 
                                   sEntry.c1, sEntry.c2, sEntry.c3);
        else
            cmap[i] = -1;
    }
  }
  else if( hBand2 == NULL && hColorMap != NULL )
  {
      int color_count;
      cmap_set = TRUE;

      color_count = MIN(256,GDALGetColorEntryCount(hColorMap));

      for(i=0; i < color_count; i++) {
          GDALColorEntry sEntry;
          
          GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );

          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] = gdTrueColorAlpha(sEntry.c1, sEntry.c2, sEntry.c3, 
                                         127 - (sEntry.c4 >> 1) );
          else
              cmap[i] = -1;
      }
  }
  else if( !truecolor && gdImg )
  {
      allocColorCube( map, gdImg, anColorCube );
  }
  else if( !truecolor )
  {
      msSetError(MS_IOERR, 
                 "Unsupported configuration for raster layer read via GDAL.",
                 "drawGDAL()");
      return -1;
  }

  /*
   * Read the entire region of interest in one gulp.  GDAL will take
   * care of downsampling and window logic.  
   */
  pabyRaw1 = (unsigned char *) malloc(dst_xsize * dst_ysize);
  if( pabyRaw1 == NULL )
  {
      msSetError(MS_MEMERR, "Allocating work image of size %dx%d failed.",
                 "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
      return -1;
  }

  if( LoadGDALImage( hBand1, 1, layer, 
                     src_xoff, src_yoff, src_xsize, src_ysize, 
                     pabyRaw1, dst_xsize, dst_ysize ) == -1 )
  {
      free( pabyRaw1 );
      return -1;
  }

  if( hBand2 != NULL && hBand3 != NULL )
  {
      pabyRaw2 = (unsigned char *) malloc(dst_xsize * dst_ysize);
      if( pabyRaw2 == NULL )
      {
          msSetError(MS_MEMERR, "Allocating work image of size %dx%d failed.",
                     "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
          return -1;
      }
      
      if( LoadGDALImage( hBand2, 2, layer, 
                         src_xoff, src_yoff, src_xsize, src_ysize, 
                         pabyRaw2, dst_xsize, dst_ysize ) == -1 )
      {
          free( pabyRaw1 );
          free( pabyRaw2 );
          return -1;
      }

      pabyRaw3 = (unsigned char *) malloc(dst_xsize * dst_ysize);
      if( pabyRaw3 == NULL )
      {
          msSetError(MS_MEMERR, "Allocating work image of size %dx%d failed.",
                     "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
          return -1;
      }
      
      if( LoadGDALImage( hBand3, 3, layer, 
                         src_xoff, src_yoff, src_xsize, src_ysize, 
                         pabyRaw3, dst_xsize, dst_ysize ) == -1 )
      {
          free( pabyRaw1 );
          free( pabyRaw2 );
          free( pabyRaw3 );
          return -1;
      }
  }

  if( hBandAlpha != NULL )
  {
      pabyRawAlpha = (unsigned char *) malloc(dst_xsize * dst_ysize);
      if( pabyRawAlpha == NULL )
      {
          msSetError(MS_MEMERR, "Allocating work image of size %dx%d failed.",
                     "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
          return -1;
      }
      
      if( LoadGDALImage( hBandAlpha, 4, layer, 
                         src_xoff, src_yoff, src_xsize, src_ysize, 
                         pabyRawAlpha, dst_xsize, dst_ysize ) == -1 )
      {
          free( pabyRaw1 );
          if( pabyRaw3 != NULL )
          {
              free( pabyRaw2 );
              free( pabyRaw3 );
          }

          free( pabyRawAlpha );
          return -1;
      }
  }
  else
      pabyRawAlpha = NULL;

/* -------------------------------------------------------------------- */
/*      Single band plus colormap with alpha blending to 8bit.          */
/* -------------------------------------------------------------------- */
  if( hBand2 == NULL && !truecolor && gdImg && hBandAlpha != NULL )
  {
      assert( cmap_set );
      k = 0;

      for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
      {
          int   result, alpha;
          
          for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ )
          {
              alpha = pabyRawAlpha[k];

              result = cmap[pabyRaw1[k++]];

              /* 
              ** We don't do alpha blending in non-truecolor mode, just
              ** threshold the point on/off at alpha=128.
              */

              if( result != -1 && alpha >= 128 )
                  gdImg->pixels[i][j] = result;
          }
      }

      assert( k == dst_xsize * dst_ysize );
  }

/* -------------------------------------------------------------------- */
/*      Single band plus colormap (no alpha) to 8bit.                   */
/* -------------------------------------------------------------------- */
  else if( hBand2 == NULL && !truecolor && gdImg )
  {
      assert( cmap_set );
      k = 0;

      for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
      {
          int   result;
          
          for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ )
          {
              result = cmap[pabyRaw1[k++]];
              if( result != -1 )
              {
                  gdImg->pixels[i][j] = result;
              }
          }
      }

      assert( k == dst_xsize * dst_ysize );
  }

/* -------------------------------------------------------------------- */
/*      Single band plus colormap and alpha to truecolor.               */
/* -------------------------------------------------------------------- */
  else if( hBand2 == NULL && truecolor && gdImg && hBandAlpha != NULL )
  {
      assert( cmap_set );

      k = 0;
      for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
      {
          int   result, alpha;
          
          for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ )
          {
              alpha = pabyRawAlpha[k];
              result = cmap[pabyRaw1[k++]];

              if( result == -1 || alpha < 2 )
                  /* do nothing - transparent */;
              else if( alpha > 253 )
                  gdImg->tpixels[i][j] = result;
              else
              {
                  /* mix alpha into "result" */
                  result += (127 - (alpha >> 1)) << 24;
                  gdImg->tpixels[i][j] = 
                      msAlphaBlend( gdImg->tpixels[i][j], result );
              }
          }
      }
  }

/* -------------------------------------------------------------------- */
/*      Single band plus colormap (no alpha) to truecolor               */
/* -------------------------------------------------------------------- */
  else if( hBand2 == NULL && truecolor && gdImg )
  {
      assert( cmap_set );

      k = 0;
      for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
      {
          int   result;
          
          for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ )
          {
              result = cmap[pabyRaw1[k++]];
              if( result != -1 )
                  gdImg->tpixels[i][j] = result;
          }
      }
  }

/* -------------------------------------------------------------------- */
/*      Input is 3 band RGB.  Alpha blending is mixed into the loop     */
/*      since this case is less commonly used and has lots of other     */
/*      overhead.                                                       */
/* -------------------------------------------------------------------- */
  else if( hBand3 != NULL && gdImg )
  {
      /* Dithered 24bit to 8bit conversion */
      if( !truecolor && CSLFetchBoolean( layer->processing, "DITHER", FALSE ) )
      {
#ifdef ENABLE_DITHER
          unsigned char *pabyDithered;

          pabyDithered = (unsigned char *) malloc(dst_xsize * dst_ysize);
          if( pabyDithered == NULL )
          {
              msSetError(MS_MEMERR, "Allocating work image of size %dx%d failed.",
                         "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
              return -1;
          }
          
          Dither24to8( pabyRaw1, pabyRaw2, pabyRaw3, pabyDithered,
                       dst_xsize, dst_ysize, image->format->transparent, 
                       map->imagecolor, gdImg );

          k = 0;
          for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
          {
              for( j = dst_xoff; j < dst_xoff + dst_xsize; j++, k++ )
              {
                  if( MS_VALID_COLOR( layer->offsite )
                      && pabyRaw1[k] == layer->offsite.red
                      && pabyRaw2[k] == layer->offsite.green
                      && pabyRaw3[k] == layer->offsite.blue )
                      continue;

                  if( pabyRawAlpha != NULL && pabyRawAlpha[k] == 0 )
                      continue;

                  gdImg->pixels[i][j] = pabyDithered[k];
              }
          }
          
          free( pabyDithered );
#else
          msSetError( MS_IMGERR, 
                      "DITHER not supported in this build.  Update to latest GDAL, and define the ENABLE_DITHER macro.", "drawGDAL()" );
          return -1;
#endif 
      }

      /* Color cubed 24bit to 8bit conversion. */
      else if( !truecolor )
      {
          k = 0;
          for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
          {
              for( j = dst_xoff; j < dst_xoff + dst_xsize; j++, k++ )
              {
                  int   cc_index;
                  
                  if( MS_VALID_COLOR( layer->offsite )
                      && pabyRaw1[k] == layer->offsite.red
                      && pabyRaw2[k] == layer->offsite.green
                      && pabyRaw3[k] == layer->offsite.blue )
                      continue;
                  
                  if( pabyRawAlpha != NULL && pabyRawAlpha[k] == 0 )
                      continue;

                  cc_index= RGB_INDEX(pabyRaw1[k],pabyRaw2[k],pabyRaw3[k]);
                  gdImg->pixels[i][j] = anColorCube[cc_index];
              }
          }
      }
      else if( truecolor )
      {
          k = 0;
          for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
          {