/******************************************************************************
 *
 * Project:  MapServer
 * Purpose:  Implementation of query operations on rasters. 
 * 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.
 ******************************************************************************
 *
 * $Log: maprasterquery.c,v $
 * Revision 1.23  2006/06/23 20:39:19  frank
 * fix time filter propogation to tileindex layer
 *
 * Revision 1.22  2006/06/22 17:32:25  frank
 * Use backtics version of LayerSetTimeFilter so that time
 * queries against layers with direct shapefile tileindexes
 * will work properly.
 *
 * Revision 1.21  2005/10/28 01:09:42  jani
 * MS RFC 3: Layer vtable architecture (bug 1477)
 *
 * Revision 1.20  2005/06/14 16:03:34  dan
 * Updated copyright date to 2005
 *
 * Revision 1.19  2005/04/21 15:09:28  julien
 * Bug1244: Replace USE_SHAPE_Z_M by USE_POINT_Z_M
 *
 * Revision 1.18  2005/04/15 19:32:33  julien
 * Bug 1103: Set the default tolerance value based on the layer type.
 *
 * Revision 1.17  2005/04/14 15:17:14  julien
 * Bug 1244: Remove Z and M from point by default to gain performance.
 *
 * Revision 1.16  2005/02/18 03:06:47  dan
 * Turned all C++ (//) comments into C comments (bug 1238)
 *
 * Revision 1.15  2004/11/15 20:35:02  dan
 * Added msLayerIsOpen() to all vector layer types (bug 1051)
 *
 * Revision 1.14  2004/11/03 16:40:24  frank
 * modified raster queries to properly set the classindex in the resultcache.
 *
 * Revision 1.13  2004/11/02 02:08:30  frank
 * Changed to use "value_list" instead of "values" for a list of all
 * pixel values.  [values] has a special substitution rule in
 * maptemplate.c, so we need to avoid it.
 *
 * Revision 1.12  2004/11/02 01:33:58  frank
 * Added logic to reproject location back to map coordiantes before doing
 * the distance or "in shape" tests.  This also ensures results are returned
 * in map coordianates.
 * http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=1021
 *
 * Revision 1.11  2004/11/01 21:32:01  frank
 * Added pointer to raster query classification bug.
 *
 * Revision 1.10  2004/10/21 04:30:56  frank
 * Added standardized headers.  Added MS_CVSID().
 *
 * Revision 1.9  2004/10/13 20:09:24  frank
 * generate a respectible error in msRASTERLayerOpen()
 *
 * Revision 1.8  2004/09/28 20:44:52  frank
 * avoid casting warning
 *
 * Revision 1.7  2004/05/18 16:52:48  frank
 * ensure things build clean without USE_GDAL
 *
 * Revision 1.6  2004/05/13 03:27:07  frank
 * fixed bug in query rect, return point geometries
 *
 * Revision 1.5  2004/05/12 20:55:49  frank
 * don't wipe query results cache between tiles
 *
 * Revision 1.4  2004/04/27 16:53:52  frank
 * added max result support
 *
 * Revision 1.3  2004/04/08 20:32:18  frank
 * now substantially complete
 *
 * Revision 1.2  2004/04/06 18:45:56  frank
 * initial implementation working
 *
 * Revision 1.1  2004/03/09 21:11:08  frank
 * New
 *
 */

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

MS_CVSID("$Id: maprasterquery.c,v 1.23 2006/06/23 20:39:19 frank Exp $")

int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, int tile, 
                          long record);

#ifdef USE_GDAL

/* ==================================================================== */
/*      For now the rasterLayerInfo lives here since it is just used    */
/*      to hold information related to queries.                         */
/* ==================================================================== */
typedef struct {

    /* query cache results */
    int query_results;
    int query_alloc_max;
    int query_request_max;
    int query_result_hard_max;
    int raster_query_mode;
    int band_count;

    int refcount;

    rectObj which_rect;
    int     next_shape;

    double *qc_x;
    double *qc_y;
    float *qc_values;
    int    *qc_class;
    int    *qc_red;
    int    *qc_green;
    int    *qc_blue;
    int    *qc_count;
    int    *qc_tileindex;

    /* query bound in force */
    shapeObj *searchshape;

    /* Only nearest result to this point. */
    int      range_mode; /* MS_SINGLE, MS_MULTI or -1 (skip test) */
    double   range_dist;
    pointObj target_point;

    GDALColorTableH hCT;

    int      current_tile;

} rasterLayerInfo;

#define RQM_UNKNOWN               0
#define RQM_ENTRY_PER_PIXEL       1
#define RQM_HIST_ON_CLASS         2
#define RQM_HIST_ON_VALUE         3

extern int InvGeoTransform( double *gt_in, double *gt_out );
#define GEO_TRANS(tr,x,y)  ((tr)[0]+(tr)[1]*(x)+(tr)[2]*(y))

/************************************************************************/
/*                             addResult()                              */
/*                                                                      */
/*      this is a copy of the code in mapquery.c.  Should we prepare    */
/*      a small public API for managing results caches?                 */
/************************************************************************/

static int addResult(resultCacheObj *cache, int classindex, int shapeindex, int tileindex)
{
  int i;

  if(cache->numresults == cache->cachesize) { /* just add it to the end */
    if(cache->cachesize == 0)
      cache->results = (resultCacheMemberObj *) malloc(sizeof(resultCacheMemberObj)*MS_RESULTCACHEINCREMENT);
    else
      cache->results = (resultCacheMemberObj *) realloc(cache->results, sizeof(resultCacheMemberObj)*(cache->cachesize+MS_RESULTCACHEINCREMENT));
    if(!cache->results) {
      msSetError(MS_MEMERR, "Realloc() error.", "addResult()");
      return(MS_FAILURE);
    }   
    cache->cachesize += MS_RESULTCACHEINCREMENT;
  }

  i = cache->numresults;

  cache->results[i].classindex = classindex;
  cache->results[i].tileindex = tileindex;
  cache->results[i].shapeindex = shapeindex;
  cache->numresults++;

  return(MS_SUCCESS);
}

/************************************************************************/
/*                       msRasterLayerInfoFree()                        */
/************************************************************************/

static void msRasterLayerInfoFree( layerObj *layer )

{
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    if( rlinfo == NULL )
        return;

    if( rlinfo->qc_x != NULL )
    {
        free( rlinfo->qc_x );
        free( rlinfo->qc_y );
    }

    if( rlinfo->qc_values )
        free( rlinfo->qc_values );

    if( rlinfo->qc_class )
    {
        free( rlinfo->qc_class );
    }

    if( rlinfo->qc_red )
    {
        free( rlinfo->qc_red );
        free( rlinfo->qc_green );
        free( rlinfo->qc_blue );
    }
    
    if( rlinfo->qc_count != NULL )
        free( rlinfo->qc_count );

    if( rlinfo->qc_tileindex != NULL )
        free( rlinfo->qc_tileindex );

    free( rlinfo );

    layer->layerinfo = NULL;
}

/************************************************************************/
/*                    msRasterLayerInfoInitialize()                     */
/************************************************************************/

static void msRasterLayerInfoInitialize( layerObj *layer )

{
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    if( rlinfo != NULL )
        return;

    rlinfo = (rasterLayerInfo *) calloc(1,sizeof(rasterLayerInfo));
    layer->layerinfo = rlinfo;
    
    rlinfo->band_count = -1;
    rlinfo->raster_query_mode = RQM_ENTRY_PER_PIXEL;
    rlinfo->range_mode = -1; /* inactive */
    rlinfo->refcount = 0;
    
    /* We need to do this or the layer->layerinfo will be interpreted */
    /* as shapefile access info because the default connectiontype is */
    /* MS_SHAPEFILE. */
    layer->connectiontype = MS_RASTER;

    rlinfo->query_result_hard_max = 1000000;

    if( CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" ) 
        != NULL )
    {
        rlinfo->query_result_hard_max = 
            atoi(CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" ));
    }
}

/************************************************************************/
/*                       msRasterQueryAddPixel()                        */
/************************************************************************/

static void msRasterQueryAddPixel( layerObj *layer, pointObj *location, 
                                   float *values )

{
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
    int red = 0, green = 0, blue = 0, nodata = FALSE;
    int p_class = -1;

    if( rlinfo->query_results == rlinfo->query_result_hard_max )
        return;

/* -------------------------------------------------------------------- */
/*      Is this our first time in?  If so, do an initial allocation     */
/*      for the data arrays suitable to our purposes.                   */
/* -------------------------------------------------------------------- */
    if( rlinfo->query_alloc_max == 0 )
    {
        rlinfo->query_alloc_max = 2;

        switch( rlinfo->raster_query_mode )
        {
          case RQM_ENTRY_PER_PIXEL:
            rlinfo->qc_x = (double *) 
                calloc(sizeof(double),rlinfo->query_alloc_max);
            rlinfo->qc_y = (double *) 
                calloc(sizeof(double),rlinfo->query_alloc_max);
            rlinfo->qc_values = (float *) 
                calloc(sizeof(float),
                       rlinfo->query_alloc_max*rlinfo->band_count);
            rlinfo->qc_red = (int *) 
                calloc(sizeof(int),rlinfo->query_alloc_max);
            rlinfo->qc_green = (int *) 
                calloc(sizeof(int),rlinfo->query_alloc_max);
            rlinfo->qc_blue = (int *) 
                calloc(sizeof(int),rlinfo->query_alloc_max);
            if( layer->numclasses > 0 )
                rlinfo->qc_class = (int *) 
                    calloc(sizeof(int),rlinfo->query_alloc_max);
            break;
           
          case RQM_HIST_ON_CLASS:
            break;
            
          case RQM_HIST_ON_VALUE:
            break;

          default:
            assert( FALSE );
        }
    }

/* -------------------------------------------------------------------- */
/*      Reallocate the data arrays larger if they are near the max      */
/*      now.                                                            */
/* -------------------------------------------------------------------- */
    if( rlinfo->query_results == rlinfo->query_alloc_max )
    {
        rlinfo->query_alloc_max = rlinfo->query_alloc_max * 2 + 100;

        if( rlinfo->qc_x != NULL )
            rlinfo->qc_x = realloc(rlinfo->qc_x,
                                   sizeof(double) * rlinfo->query_alloc_max);
        if( rlinfo->qc_y != NULL )
            rlinfo->qc_y = realloc(rlinfo->qc_y,
                                   sizeof(double) * rlinfo->query_alloc_max);
        if( rlinfo->qc_values != NULL )
            rlinfo->qc_values = 
                realloc(rlinfo->qc_values,
                        sizeof(float) * rlinfo->query_alloc_max 
                        * rlinfo->band_count );
        if( rlinfo->qc_class != NULL )
            rlinfo->qc_class = realloc(rlinfo->qc_class,
                                       sizeof(int) * rlinfo->query_alloc_max);
        if( rlinfo->qc_red != NULL )
            rlinfo->qc_red = realloc(rlinfo->qc_red,
                                       sizeof(int) * rlinfo->query_alloc_max);
        if( rlinfo->qc_green != NULL )
            rlinfo->qc_green = realloc(rlinfo->qc_green,
                                       sizeof(int) * rlinfo->query_alloc_max);
        if( rlinfo->qc_blue != NULL )
            rlinfo->qc_blue = realloc(rlinfo->qc_blue,
                                       sizeof(int) * rlinfo->query_alloc_max);
        if( rlinfo->qc_count != NULL )
            rlinfo->qc_count = realloc(rlinfo->qc_count,
                                       sizeof(int) * rlinfo->query_alloc_max);
        if( rlinfo->qc_count != NULL )
            rlinfo->qc_tileindex = realloc(rlinfo->qc_tileindex,
                                       sizeof(int) * rlinfo->query_alloc_max);
    }

/* -------------------------------------------------------------------- */
/*      Handle classification.                                          */
/*                                                                      */
/*      NOTE: The following is really quite inadequate to deal with     */
/*      classifications based on [red], [green] and [blue] as           */
/*      described in:                                                   */
/*       http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=1021         */
/* -------------------------------------------------------------------- */
    if( rlinfo->qc_class != NULL )
    {
        p_class = msGetClass_Float(layer, values[0] );

        if( p_class == -1 )
            nodata = TRUE;
        else
        {
            rlinfo->qc_class[rlinfo->query_results] = p_class;
            red   = layer->class[p_class].styles[0].color.red;
            green = layer->class[p_class].styles[0].color.green;
            blue  = layer->class[p_class].styles[0].color.blue;
        }
    }

/* -------------------------------------------------------------------- */
/*      Handle colormap                                                 */
/* -------------------------------------------------------------------- */
    else if( rlinfo->hCT != NULL )
    {
        int pct_index = (int) floor(values[0]);
        GDALColorEntry sEntry;

        if( GDALGetColorEntryAsRGB( rlinfo->hCT, pct_index, &sEntry ) )
        {
            red = sEntry.c1;
            green = sEntry.c2;
            blue = sEntry.c3;

            if( sEntry.c4 == 0 )
                nodata = TRUE;
        }
        else
            nodata = TRUE;
    }

/* -------------------------------------------------------------------- */
/*      Color derived from greyscale value.                             */
/* -------------------------------------------------------------------- */
    else
    {
        if( rlinfo->band_count >= 3 )
        {
            red = (int) MAX(0,MIN(255,values[0]));
            green = (int) MAX(0,MIN(255,values[1]));
            blue = (int) MAX(0,MIN(255,values[2]));
        }
        else
        {
            red = green = blue = (int) MAX(0,MIN(255,values[0]));
        }
    }

/* -------------------------------------------------------------------- */
/*      Record the color.                                               */
/* -------------------------------------------------------------------- */
    rlinfo->qc_red[rlinfo->query_results] = red;
    rlinfo->qc_green[rlinfo->query_results] = green;
    rlinfo->qc_blue[rlinfo->query_results] = blue;

/* -------------------------------------------------------------------- */
/*      Record spatial location.                                        */
/* -------------------------------------------------------------------- */
    if( rlinfo->qc_x != NULL )
    {
        rlinfo->qc_x[rlinfo->query_results] = location->x;
        rlinfo->qc_y[rlinfo->query_results] = location->y;
    }

/* -------------------------------------------------------------------- */
/*      Record actual pixel value(s).                                   */
/* -------------------------------------------------------------------- */
    if( rlinfo->qc_values != NULL )
        memcpy( rlinfo->qc_values + rlinfo->query_results * rlinfo->band_count,
                values, sizeof(float) * rlinfo->band_count );

/* -------------------------------------------------------------------- */
/*      Add to the results cache.                                       */
/* -------------------------------------------------------------------- */
    if( ! nodata )
    {
        addResult( layer->resultcache, p_class, rlinfo->query_results, 0 );
        rlinfo->query_results++;
    }
}

/************************************************************************/
/*                       msRasterQueryByRectLow()                       */
/************************************************************************/

static int 
msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
                       rectObj queryRect) 

{
    double	adfGeoTransform[6], adfInvGeoTransform[6];
    double      dfXMin, dfYMin, dfXMax, dfYMax, dfX, dfY, dfAdjustedRange;
    int         nWinXOff, nWinYOff, nWinXSize, nWinYSize;
    int         nRXSize, nRYSize;
    float       *pafRaster;
    int         nBandCount, *panBandMap, iPixel, iLine;
    CPLErr      eErr;
    rasterLayerInfo *rlinfo;
    rectObj     searchrect;
    int         needReproject = MS_FALSE;

    rlinfo = (rasterLayerInfo *) layer->layerinfo;

/* -------------------------------------------------------------------- */
/*      Reproject the search rect into the projection of this           */
/*      layer/file if needed.                                           */
/* -------------------------------------------------------------------- */
    searchrect = queryRect;
#ifdef USE_PROJ
    if(layer->project 
       && msProjectionsDiffer(&(layer->projection), &(map->projection)))
    {
        msProjectRect(&(map->projection), &(layer->projection), &searchrect);
        needReproject = MS_TRUE;
    }
    else
        layer->project = MS_FALSE;
#endif

/* -------------------------------------------------------------------- */
/*      Transform the rectangle in target ground coordinates to         */
/*      pixel/line extents on the file.  Process all 4 corners, to      */
/*      build extents.                                                  */
/* -------------------------------------------------------------------- */
    nRXSize = GDALGetRasterXSize( hDS );
    nRYSize = GDALGetRasterYSize( hDS );

    msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
    InvGeoTransform( adfGeoTransform, adfInvGeoTransform );

    /* top left */
    dfXMin = dfXMax = GEO_TRANS(adfInvGeoTransform,
                                searchrect.minx, searchrect.maxy);
    dfYMin = dfYMax = GEO_TRANS(adfInvGeoTransform+3,
                                searchrect.minx, searchrect.maxy);
                                
    /* top right */
    dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.maxx, searchrect.maxy);
    dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.maxy);
    dfXMin = MIN(dfXMin,dfX);
    dfXMax = MAX(dfXMax,dfX);
    dfYMin = MIN(dfYMin,dfY);
    dfYMax = MAX(dfYMax,dfY);
    
    /* bottom left */
    dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.minx, searchrect.miny);
    dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.minx, searchrect.miny);
    dfXMin = MIN(dfXMin,dfX);
    dfXMax = MAX(dfXMax,dfX);
    dfYMin = MIN(dfYMin,dfY);
    dfYMax = MAX(dfYMax,dfY);
    
    /* bottom right */
    dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.maxx, searchrect.miny);
    dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.miny);
    dfXMin = MIN(dfXMin,dfX);
    dfXMax = MAX(dfXMax,dfX);
    dfYMin = MIN(dfYMin,dfY);
    dfYMax = MAX(dfYMax,dfY);

/* -------------------------------------------------------------------- */
/*      Trim the rectangle to the area of the file itself, but out      */
/*      to the edges of the touched edge pixels.                        */
/* -------------------------------------------------------------------- */
    dfXMin = MAX(0.0,floor(dfXMin));
    dfYMin = MAX(0.0,floor(dfYMin));
    dfXMax = MIN(nRXSize,ceil(dfXMax));
    dfYMax = MIN(nRYSize,ceil(dfYMax));

/* -------------------------------------------------------------------- */
/*      Convert to integer offset/size values.                          */
/* -------------------------------------------------------------------- */
    nWinXOff = (int) dfXMin;
    nWinYOff = (int) dfYMin;
    nWinXSize = (int) (dfXMax - dfXMin);
    nWinYSize = (int) (dfYMax - dfYMin);

/* -------------------------------------------------------------------- */
/*      What bands are we operating on?                                 */
/* -------------------------------------------------------------------- */
    panBandMap = msGetGDALBandList( layer, hDS, 0, &nBandCount );
    
    if( rlinfo->band_count == -1 )
        rlinfo->band_count = nBandCount;

    if( nBandCount != rlinfo->band_count )
    {
        msSetError( MS_IMGERR, 
                    "Got %d bands, but expected %d bands.", 
                    "msRasterQueryByRectLow()", 
                    nBandCount, rlinfo->band_count );

        return -1;
    }

/* -------------------------------------------------------------------- */
/*      Try to load the raster data.  For now we just load the first    */
/*      band in the file.  Later we will deal with the various band     */
/*      selection criteria.                                             */
/* -------------------------------------------------------------------- */
    pafRaster = (float *) 
        calloc(sizeof(float),nWinXSize*nWinYSize*nBandCount);

#if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1199
    eErr = GDALDatasetRasterIO( hDS, GF_Read, 
                                nWinXOff, nWinYOff, nWinXSize, nWinYSize,
                                pafRaster, nWinXSize, nWinYSize, GDT_Float32,
                                nBandCount, panBandMap, 0, 0, 0 );
#else
    /*
     * The above could actually be implemented for pre-1.2.0 GDALs
     * reading band by band, but it would be hard to do and test and would
     * be very rarely useful so we skip it.
     */
    msSetError(MS_IMGERR, 
               "Raster query support requires GDAL 1.2.0 or newer.", 
               "msRasterQueryByRectLow()" );
    free( pafRaster );
    return MS_FAILURE;
#endif

    
    if( eErr != CE_None )
    {
        msSetError( MS_IOERR, "GDALDatasetRasterIO() failed: %s", 
                    CPLGetLastErrorMsg(), 
                    "msRasterQueryByRectLow()" );

        free( pafRaster );
        return -1;
    }

/* -------------------------------------------------------------------- */
/*      Fetch color table for intepreting colors if needed.             */
/* -------------------------------------------------------------------- */
    rlinfo->hCT = GDALGetRasterColorTable( 
        GDALGetRasterBand( hDS, panBandMap[0] ) );

    free( panBandMap );

/* -------------------------------------------------------------------- */
/*      When computing whether pixels are within range we do it         */
/*      based on the center of the pixel to the target point but        */
/*      really it ought to be the nearest point on the pixel.  It       */
/*      would be too much trouble to do this rigerously, so we just     */
/*      add a fudge factor so that a range of zero will find the        */
/*      pixel the target falls in at least.                             */
/* -------------------------------------------------------------------- */
    dfAdjustedRange = 
        sqrt(adfGeoTransform[1] * adfGeoTransform[1]
             + adfGeoTransform[2] * adfGeoTransform[2]
             + adfGeoTransform[4] * adfGeoTransform[4]
             + adfGeoTransform[5] * adfGeoTransform[5]) * 0.5 * 1.41421356237
        + sqrt( rlinfo->range_dist );
    dfAdjustedRange = dfAdjustedRange * dfAdjustedRange;

/* -------------------------------------------------------------------- */
/*      Loop over all pixels determining which are "in".                */
/* -------------------------------------------------------------------- */
    for( iLine = 0; iLine < nWinYSize; iLine++ )
    {
        for( iPixel = 0; iPixel < nWinXSize; iPixel++ )
        {
            pointObj  sPixelLocation;

            if( rlinfo->query_results == rlinfo->query_result_hard_max )
                break;

            /* transform pixel/line to georeferenced */
            sPixelLocation.x = 
                GEO_TRANS(adfGeoTransform, 
                          iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );
            sPixelLocation.y = 
                GEO_TRANS(adfGeoTransform+3, 
                          iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );

            /* If projections differ, convert this back into the map  */
            /* projection for distance testing, and comprison to the  */
            /* search shape.  */
            if( needReproject )
                msProjectPoint( &(layer->projection), &(map->projection), 
                                &sPixelLocation );

            if( rlinfo->searchshape != NULL
                && msIntersectPointPolygon( &sPixelLocation, 
                                            rlinfo->searchshape ) == MS_FALSE )
                continue;

            if( rlinfo->range_mode >= 0 )
            {
                double dist;

                dist = (rlinfo->target_point.x - sPixelLocation.x) 
                    * (rlinfo->target_point.x - sPixelLocation.x) 
                    + (rlinfo->target_point.y - sPixelLocation.y)
                    * (rlinfo->target_point.y - sPixelLocation.y);

                if( dist >= dfAdjustedRange )
                    continue;

                /* If we can only have one feature, trim range and clear */
                /* previous result.  */
                if( rlinfo->range_mode == MS_SINGLE )
                {
                    rlinfo->range_dist = dist;
                    rlinfo->query_results = 0;
                }
            }

            msRasterQueryAddPixel( layer, &sPixelLocation, 
                                   pafRaster 
                                   + (iLine*nWinXSize + iPixel) * nBandCount );
        }
    }

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
/* -------------------------------------------------------------------- */
    free( pafRaster );

    return MS_SUCCESS;
}
#endif /* def USE_GDAL */

/************************************************************************/
/*                        msRasterQueryByRect()                         */
/************************************************************************/

int msRasterQueryByRect(mapObj *map, layerObj *layer, rectObj queryRect) 

{
#ifndef USE_GDAL
    msSetError( MS_IMGERR, 
                "Rasters queries only supported with GDAL support enabled.",
                "msRasterQueryByRect()" );
    return MS_FAILURE;
#else
    int status = MS_SUCCESS;
    char *filename=NULL;

    int t;
    int tileitemindex=-1;
    shapefileObj tilefile;
    char tilename[MS_PATH_LENGTH];
    int numtiles=1; /* always at least one tile */
    char szPath[MS_MAXPATHLEN];
    rectObj searchrect;
    rasterLayerInfo *rlinfo = NULL;

/* -------------------------------------------------------------------- */
/*      Get the layer info.                                             */
/* -------------------------------------------------------------------- */
    msRasterLayerInfoInitialize( layer );
    rlinfo = (rasterLayerInfo *) layer->layerinfo;

/* -------------------------------------------------------------------- */
/*      Clear old results cache.                                        */
/* -------------------------------------------------------------------- */
    if(layer->resultcache) {
      if(layer->resultcache->results) free(layer->resultcache->results);
      free(layer->resultcache);
      layer->resultcache = NULL;
    }

/* -------------------------------------------------------------------- */
/*      Initialize the results cache.                                   */
/* -------------------------------------------------------------------- */
    layer->resultcache = (resultCacheObj *)malloc(sizeof(resultCacheObj)); 
    layer->resultcache->results = NULL;
    layer->resultcache->numresults = layer->resultcache->cachesize = 0;
    layer->resultcache->bounds.minx =
        layer->resultcache->bounds.miny =
        layer->resultcache->bounds.maxx =
        layer->resultcache->bounds.maxy = -1;

/* -------------------------------------------------------------------- */
/*      Check if we should really be acting on this layer and           */
/*      provide debug info in various cases.                            */
/* -------------------------------------------------------------------- */
    if(layer->debug > 0 || map->debug > 1)
        msDebug( "msRasterQueryByRect(%s): entering.\n", layer->name );

    if(!layer->data && !layer->tileindex) {
        if(layer->debug > 0 || map->debug > 0 )
            msDebug( "msRasterQueryByRect(%s): layer data and tileindex NULL ... doing nothing.", layer->name );
        return MS_SUCCESS;
    }

    if((layer->status != MS_ON) && layer->status != MS_DEFAULT ) {
        if(layer->debug > 0 )
            msDebug( "msRasterQueryByRect(%s): not status ON or DEFAULT, doing nothing.", layer->name );
        return MS_SUCCESS;
    }

/* -------------------------------------------------------------------- */
/*      Open tile index if we have one.                                 */
/* -------------------------------------------------------------------- */
    if(layer->tileindex) { /* we have in index file */
        if(msSHPOpenFile(&tilefile, "rb", 
                         msBuildPath3(szPath, map->mappath, map->shapepath, 
                                      layer->tileindex)) == -1) 
            if(msSHPOpenFile(&tilefile, "rb", msBuildPath(szPath, map->mappath, layer->tileindex)) == -1) 
                return(MS_FAILURE);    

        tileitemindex = msDBFGetItemIndex(tilefile.hDBF, layer->tileitem);
        if(tileitemindex == -1) 
            return(MS_FAILURE);

        searchrect = queryRect;

#ifdef USE_PROJ
        if((map->projection.numargs > 0) && (layer->projection.numargs > 0))
            msProjectRect(&map->projection, &layer->projection, &searchrect); /* project the searchrect to source coords */
#endif
        status = msSHPWhichShapes(&tilefile, searchrect, layer->debug);
        if(status != MS_SUCCESS) 
            numtiles = 0; /* could be MS_DONE or MS_FAILURE */
        else
            numtiles = tilefile.numshapes;
    }

/* -------------------------------------------------------------------- */
/*      Iterate over all tiles (just one in untiled case).              */
/* -------------------------------------------------------------------- */
    for(t=0; t<numtiles && status == MS_SUCCESS; t++) { 

        GDALDatasetH  hDS;

        rlinfo->current_tile = t;

/* -------------------------------------------------------------------- */
/*      Get filename.                                                   */
/* -------------------------------------------------------------------- */
        if(layer->tileindex) {
            if(!msGetBit(tilefile.status,t)) continue; /* on to next tile */
            if(layer->data == NULL) /* assume whole filename is in attribute field */
                filename = (char*)msDBFReadStringAttribute(tilefile.hDBF, t, tileitemindex);
            else {  
                sprintf(tilename,"%s/%s", msDBFReadStringAttribute(tilefile.hDBF, t, tileitemindex) , layer->data);
                filename = tilename;
            }
        } else {
            filename = layer->data;
        }

        if(strlen(filename) == 0) continue;

/* -------------------------------------------------------------------- */
/*      Open the file.                                                  */
/* -------------------------------------------------------------------- */
        msGDALInitialize();

        msAcquireLock( TLOCK_GDAL );
        hDS = GDALOpen(
            msBuildPath3( szPath, map->mappath, map->shapepath, filename ), 
            GA_ReadOnly );
        
        if( hDS == NULL )
        {
            msReleaseLock( TLOCK_GDAL );

#ifndef IGNORE_MISSING_DATA
            if( layer->debug || map->debug )
                msSetError( MS_IMGERR, 
                            "Unable to open file %s for layer %s ... fatal error.\n%s", 
                            szPath, layer->name, CPLGetLastErrorMsg(),
                            "msRasterQueryByRect()" );
            
            return(MS_FAILURE);
#else
            if( layer->debug || map->debug )
                msDebug( "Unable to open file %s for layer %s ... ignoring this missing data.\n%s", 
                         filename, layer->name, CPLGetLastErrorMsg() );
            
#endif
            continue;
        }

/* -------------------------------------------------------------------- */
/*      Update projectionObj if AUTO.                                   */
/* -------------------------------------------------------------------- */
        if (layer->projection.numargs > 0 && 
            EQUAL(layer->projection.args[0], "auto"))
        {
            const char *pszWKT;

            pszWKT = GDALGetProjectionRef( hDS );

            if( pszWKT != NULL && strlen(pszWKT) > 0 )
            {
                if( msOGCWKT2ProjectionObj(pszWKT, &(layer->projection),
                                           layer->debug ) != MS_SUCCESS )
                {
                    char	szLongMsg[MESSAGELENGTH*2];
                    errorObj *ms_error = msGetErrorObj();

                    sprintf( szLongMsg, 
                             "%s\n"
                             "PROJECTION AUTO cannot be used for this "
                             "GDAL raster (`%s').",
                             ms_error->message, filename);
                    szLongMsg[MESSAGELENGTH-1] = '\0';

                    msSetError(MS_OGRERR, "%s","msDrawRasterLayer()",
                               szLongMsg);

                    msReleaseLock( TLOCK_GDAL );
                    return(MS_FAILURE);
                }
            }
        }

/* -------------------------------------------------------------------- */
/*      Perform actual query on this file.                              */
/* -------------------------------------------------------------------- */
        if( status == MS_SUCCESS )
            status = msRasterQueryByRectLow( map, layer, hDS, queryRect );

        GDALClose( hDS );
        msReleaseLock( TLOCK_GDAL );

    } /* next tile */

    if(layer->tileindex) /* tiling clean-up */
        msSHPCloseFile(&tilefile);    

    return status;
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                        msRasterQueryByShape()                        */
/************************************************************************/

int msRasterQueryByShape(mapObj *map, layerObj *layer, shapeObj *selectshape)

{
#ifndef USE_GDAL
    msSetError( MS_IMGERR, 
                "Rasters queries only supported with GDAL support enabled.",
                "msRasterQueryByRect()" );
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = NULL;
    int status;

    msRasterLayerInfoInitialize( layer );

    rlinfo = (rasterLayerInfo *) layer->layerinfo;
    rlinfo->searchshape = selectshape;

    msComputeBounds( selectshape );
    status = msRasterQueryByRect( map, layer, selectshape->bounds );

    rlinfo->searchshape = NULL;

    return status;
#endif /* USE_GDAL  */
}

/************************************************************************/
/*                        msRasterQueryByPoint()                        */
/************************************************************************/

int msRasterQueryByPoint(mapObj *map, layerObj *layer, int mode, pointObj p, 
                         double buffer)
{
#ifndef USE_GDAL
    msSetError( MS_IMGERR, 
                "Rasters queries only supported with GDAL support enabled.",
                "msRasterQueryByRect()" );
    return MS_FAILURE;
#else
    int result;
    double layer_tolerance;
    rectObj bufferRect;
    rasterLayerInfo *rlinfo = NULL;

    msRasterLayerInfoInitialize( layer );

    rlinfo = (rasterLayerInfo *) layer->layerinfo;

/* -------------------------------------------------------------------- */
/*      If the buffer is not set, then use layer tolerances             */
/*      instead.   The "buffer" distince is now in georeferenced        */
/*      units.  Note that tolerances in pixels are basically map        */
/*      display pixels, not underlying raster pixels.  It isn't         */
/*      clear that there is any way of requesting a buffer size in      */
/*      underlying pixels.                                              */
/* -------------------------------------------------------------------- */
    if(buffer <= 0) { /* use layer tolerance */

        /* Get the layer tolerance
           default is 3 for point and line layers, 0 for others */
        if(layer->tolerance == -1)
            if(layer->status == MS_LAYER_POINT || 
               layer->status == MS_LAYER_LINE)
                layer_tolerance = 3;
            else
                layer_tolerance = 0;
        else
            layer_tolerance = layer->tolerance;

        if(layer->toleranceunits == MS_PIXELS)
            buffer = layer_tolerance 
                * msAdjustExtent(&(map->extent), map->width, map->height);
        else
            buffer = layer_tolerance 
                * (msInchesPerUnit(layer->toleranceunits,0)
                   / msInchesPerUnit(map->units,0));
    }

/* -------------------------------------------------------------------- */
/*      Setup target point information, at this point they are in       */
/*      map coordinates.                                                */
/* -------------------------------------------------------------------- */
    rlinfo->range_dist = buffer * buffer;
    rlinfo->target_point = p;

/* -------------------------------------------------------------------- */
/*      if we are in the MS_SINGLE mode, first try a query with zero    */
/*      tolerance.  If this gets a raster pixel then we can be          */
/*      reasonably assured that it is the closest to the query          */
/*      point.  This will potentially be must more efficient than       */
/*      processing all pixels within the tolerance.                     */
/* -------------------------------------------------------------------- */
    if( mode == MS_SINGLE )
    {
        rectObj pointRect;

        pointRect.minx = p.x;
        pointRect.maxx = p.x;
        pointRect.miny = p.y;
        pointRect.maxy = p.y;

        rlinfo->range_mode = MS_SINGLE;
        result = msRasterQueryByRect( map, layer, pointRect );

        if( rlinfo->query_results > 0 )
            return result;
    }

/* -------------------------------------------------------------------- */
/*      Setup a rectangle that is everything within the designated      */
/*      range and do a search against that.                             */
/* -------------------------------------------------------------------- */
    bufferRect.minx = p.x - buffer;
    bufferRect.maxx = p.x + buffer;
    bufferRect.miny = p.y - buffer;
    bufferRect.maxy = p.y + buffer;

    rlinfo->range_mode = mode;
    result = msRasterQueryByRect( map, layer, bufferRect );

    return result;
#endif /* USE_GDAL  */
}

/************************************************************************/
/* ==================================================================== */
/*                          RASTER Query Layer                          */
/*                                                                      */
/*      The results of a raster query are treated as a set of shapes    */
/*      that can be accessed using the normal layer semantics.          */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                         msRASTERLayerOpen()                          */
/************************************************************************/

int msRASTERLayerOpen(layerObj *layer) 
{
#ifndef USE_GDAL
    msSetError( MS_IMGERR, 
                "Rasters queries only supported with GDAL support enabled.",
                "msRasterQueryByRect()" );
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    if( rlinfo != NULL )
        rlinfo->refcount = rlinfo->refcount + 1;

    if( rlinfo == NULL )
    {
        msSetError( MS_MISCERR, 
                    "Attempt to open a RASTER layer, but this is"
                    " only supported after a raster query.",
                    "msRASTERLayerOpen()" );
        return MS_FAILURE;
    }
    else 
        return MS_SUCCESS;
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                         msRASTERIsLayerOpen()                        */
/************************************************************************/

int msRASTERLayerIsOpen(layerObj *layer) 
{
#ifndef USE_GDAL
    msSetError( MS_IMGERR, 
                "Rasters queries only supported with GDAL support enabled.",
                "msRasterLayerIsOpen()" );
    return MS_FALSE;
#else
    if (layer->layerinfo)
        return MS_TRUE;
    return MS_FALSE;
#endif
}


/************************************************************************/
/*                     msRASTERLayerFreeItemInfo()                      */
/************************************************************************/
void msRASTERLayerFreeItemInfo(layerObj *layer) 
	{}

/************************************************************************/
/*                     msRASTERLayerInitItemInfo()                      */
/*                                                                      */
/*      Perhaps we should be validating the requested items here?       */
/************************************************************************/

int msRASTERLayerInitItemInfo(layerObj *layer) 
	{ return MS_SUCCESS; }

/************************************************************************/
/*                      msRASTERLayerWhichShapes()                      */
/************************************************************************/
int msRASTERLayerWhichShapes(layerObj *layer, rectObj rect) 

{
#ifndef USE_GDAL
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    rlinfo->which_rect = rect;
    rlinfo->next_shape = 0;

    return MS_SUCCESS;
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                         msRASTERLayerClose()                         */
/************************************************************************/

int msRASTERLayerClose(layerObj *layer) 
{
#ifndef USE_GDAL
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    if( rlinfo != NULL )
    {
        rlinfo->refcount--;

        if( rlinfo->refcount < 0 )
            msRasterLayerInfoFree( layer );
    }
    return MS_SUCCESS;
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerNextShape()                       */
/************************************************************************/

int msRASTERLayerNextShape(layerObj *layer, shapeObj *shape)
{
#ifndef USE_GDAL
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    if( rlinfo->next_shape < 0 
        || rlinfo->next_shape >= rlinfo->query_results )
    {
        msFreeShape(shape);
        shape->type = MS_SHAPE_NULL;
        return MS_DONE;
    }
    else 
        return msRASTERLayerGetShape( layer, shape, 0, rlinfo->next_shape++ );
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerGetShape()                        */
/************************************************************************/

int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, int tile, 
                          long record)

{
#ifndef USE_GDAL
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
    int i;

    msFreeShape(shape);
    shape->type = MS_SHAPE_NULL;

/* -------------------------------------------------------------------- */
/*      Validate requested record id.                                   */
/* -------------------------------------------------------------------- */
    if( record < 0 || record >= rlinfo->query_results )
    {
        msSetError(MS_MISCERR, 
                   "Out of range shape index requested.  Requested %d\n"
                   "but only %d shapes available.",
                   "msRASTERLayerGetShape()",
                   record, rlinfo->query_results );
        return MS_FAILURE;
    }

/* -------------------------------------------------------------------- */
/*      Apply the geometry.                                             */
/* -------------------------------------------------------------------- */
    if( rlinfo->qc_x != NULL )
    {
	lineObj	line;
        pointObj point;

        shape->type = MS_SHAPE_POINT;
        
        line.numpoints = 1;
        line.point = &point;
        
        point.x = rlinfo->qc_x[record];
        point.y = rlinfo->qc_y[record];
#ifdef USE_POINT_Z_M
        point.m = 0.0;
#endif

        msAddLine( shape, &line );
    }

/* -------------------------------------------------------------------- */
/*      Apply the requested items.                                      */
/* -------------------------------------------------------------------- */
    if( layer->numitems > 0 )
    {
        shape->values = (char **) malloc(sizeof(char *) * layer->numitems);
        shape->numvalues = layer->numitems;
        
        for( i = 0; i < layer->numitems; i++ )
        {
            char szWork[1000];

            szWork[0] = '\0';
            if( EQUAL(layer->items[i],"x") && rlinfo->qc_x )
                sprintf( szWork, "%.8g", rlinfo->qc_x[record] );
            else if( EQUAL(layer->items[i],"y") && rlinfo->qc_y )
                sprintf( szWork, "%.8g", rlinfo->qc_y[record] );

            else if( EQUAL(layer->items[i],"value_list") && rlinfo->qc_values )
            {
                int iValue;

                for( iValue = 0; iValue < rlinfo->band_count; iValue++ )
                {
                    if( iValue != 0 )
                        strcat( szWork, "," );

                    sprintf( szWork+strlen(szWork), "%.8g", 
                             rlinfo->qc_values[record * rlinfo->band_count
                                               + iValue] );
                }
            }
            else if( EQUALN(layer->items[i],"value_",6) && rlinfo->qc_values )
            {
                int iValue = atoi(layer->items[i]+6);
                sprintf( szWork, "%.8g", 
                         rlinfo->qc_values[record*rlinfo->band_count+iValue] );
            }
            else if( EQUAL(layer->items[i],"class") && rlinfo->qc_class ) 
            {
                int p_class = rlinfo->qc_class[record];
                if( layer->class[p_class].name != NULL )
                    sprintf( szWork, "%.999s", layer->class[p_class].name );
                else
                    sprintf( szWork, "%d", p_class );
            }
            else if( EQUAL(layer->items[i],"red") && rlinfo->qc_red )
                sprintf( szWork, "%d", rlinfo->qc_red[record] );
            else if( EQUAL(layer->items[i],"green") && rlinfo->qc_green )
                sprintf( szWork, "%d", rlinfo->qc_green[record] );
            else if( EQUAL(layer->items[i],"blue") && rlinfo->qc_blue )
                sprintf( szWork, "%d", rlinfo->qc_blue[record] );
            else if( EQUAL(layer->items[i],"count") && rlinfo->qc_count )
                sprintf( szWork, "%d", rlinfo->qc_count[record] );

            shape->values[i] = strdup(szWork);
        }
    }

/* -------------------------------------------------------------------- */
/*      Eventually we should likey apply the geometry properly but      */
/*      we don't really care about the geometry for query purposes.     */
/* -------------------------------------------------------------------- */

    return MS_SUCCESS;
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerGetExtent()                       */
/************************************************************************/

int msRASTERLayerGetExtent(layerObj *layer, rectObj *extent)
	{ return MS_FAILURE; }

/************************************************************************/
/*                       msRASTERLayerGetItems()                        */
/************************************************************************/

int msRASTERLayerGetItems(layerObj *layer)
{
#ifndef USE_GDAL
    return MS_FAILURE;
#else
    rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

    layer->items = (char **) calloc(sizeof(char *),10);

    layer->numitems = 0;
    if( rlinfo->qc_x )
        layer->items[layer->numitems++] = strdup("x");
    if( rlinfo->qc_y )
        layer->items[layer->numitems++] = strdup("y");
    if( rlinfo->qc_values )
    {
        int i;
        for( i = 0; i < rlinfo->band_count; i++ )
        {
            char szName[100];
            sprintf( szName, "value_%d", i );
            layer->items[layer->numitems++] = strdup(szName);
        }
        layer->items[layer->numitems++] = strdup("value_list");
    }
    if( rlinfo->qc_class )
        layer->items[layer->numitems++] = strdup("class");
    if( rlinfo->qc_red )
        layer->items[layer->numitems++] = strdup("red");
    if( rlinfo->qc_green )
        layer->items[layer->numitems++] = strdup("green");
    if( rlinfo->qc_blue )
        layer->items[layer->numitems++] = strdup("blue");
    if( rlinfo->qc_count )
        layer->items[layer->numitems++] = strdup("count");

    return msRASTERLayerInitItemInfo(layer);
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                     msRASTERLayerSetTimeFilter()                     */
/*                                                                      */
/*      This function is actually just used in the context of           */
/*      setting a filter on the tileindex for time based queries.       */
/*      For instance via WMS requests.  So it isn't really related      */
/*      to the "raster query" support at all.                           */
/*                                                                      */
/*      If a local shapefile tileindex is in use, we will set a         */
/*      backtics filter (shapefile compatible).  If another layer is    */
/*      being used as the tileindex then we will forward the            */
/*      SetTimeFilter call to it.  If there is no tileindex in          */
/*      place, we do nothing.                                           */
/************************************************************************/

int msRASTERLayerSetTimeFilter(layerObj *layer, const char *timestring, 
                               const char *timefield)
{
    int tilelayerindex;

/* -------------------------------------------------------------------- */
/*      If we don't have a tileindex the time filter has no effect.     */
/* -------------------------------------------------------------------- */
    if( layer->tileindex == NULL ) 
        return MS_SUCCESS;

/* -------------------------------------------------------------------- */
/*      Find the tileindex layer.                                       */
/* -------------------------------------------------------------------- */
    tilelayerindex = msGetLayerIndex(layer->map, layer->tileindex);

/* -------------------------------------------------------------------- */
/*      If we are using a local shapefile as our tileindex (that is     */
/*      to say, the tileindex name is not of another layer), then we    */
/*      just install a backtics style filter on the raster layer.       */
/*      This is propogated to the "working layer" created for the       */
/*      tileindex by code in mapraster.c.                               */
/* -------------------------------------------------------------------- */
    if( tilelayerindex == -1 )
        return msLayerMakeBackticsTimeFilter( layer, timestring, timefield );

/* -------------------------------------------------------------------- */
/*      Otherwise we invoke the tileindex layers SetTimeFilter          */
/*      method.                                                         */
/* -------------------------------------------------------------------- */
    return msLayerSetTimeFilter( layer->map->layers + tilelayerindex,
                                 timestring, timefield );
}

/************************************************************************/
/*                msRASTERLayerInitializeVirtualTable()                 */
/************************************************************************/

int
msRASTERLayerInitializeVirtualTable(layerObj *layer)
{
    assert(layer != NULL);
    assert(layer->vtable != NULL);

    layer->vtable->LayerInitItemInfo = msRASTERLayerInitItemInfo;
    layer->vtable->LayerFreeItemInfo = msRASTERLayerFreeItemInfo;
    layer->vtable->LayerOpen = msRASTERLayerOpen;
    layer->vtable->LayerIsOpen = msRASTERLayerIsOpen;
    layer->vtable->LayerWhichShapes = msRASTERLayerWhichShapes;
    layer->vtable->LayerNextShape = msRASTERLayerNextShape;
    layer->vtable->LayerGetShape = msRASTERLayerGetShape;

    layer->vtable->LayerClose = msRASTERLayerClose;
    layer->vtable->LayerGetItems = msRASTERLayerGetItems;
    layer->vtable->LayerGetExtent = msRASTERLayerGetExtent;

    /* layer->vtable->LayerGetAutoStyle, use default */
    /* layer->vtable->LayerApplyFilterToLayer, use default */

    layer->vtable->LayerCloseConnection = msRASTERLayerClose;

    /* we use backtics for proper tileindex shapefile functioning */
    layer->vtable->LayerSetTimeFilter = msRASTERLayerSetTimeFilter;

    /* layer->vtable->LayerCreateItems, use default */
    /* layer->vtable->LayerGetNumFeatures, use default */

    return MS_SUCCESS;
}