[gdal-dev] RasterIO decimation question ?

[Frank Warmerdam]

On Thu, Mar 18, 2010 at 9:46 AM, tle <thanhvu.le at univ-pau.fr> wrote:
>
> Hi,
>
> I'm writing a terrain visualization prototype and I'm using GDAL for reading
> DEM data (GeoTif format without compressing and overview)
> When the source and destination image dimension are different, I have
> noticed a small translation of destination image to the left and bottom when
> the dimension is getting smaller and smaller.
> I have written a test program to find out what happened
>
> This is an example of how i did :
> Here is my source image with only 4x4 pixels :
> 0,0    1,0    2,0    3,0
> 0,1    1,1    2,1    3,1
> 0,2    1,2    2,2    3,2
> 0,3    1,3    2,3    3,3
>
> By performing a rasterIO with destination dimension = 4, I got :
> 0,0    1,0    2,0    3,0
> 0,1    1,1    2,1    3,1
> 0,2    1,2    2,2    3,2
> 0,3    1,3    2,3    3,3
> which is normal
>
> Now I keep the same view and reduce the destination by haft (destination
> dimension = 2), I got :
> 1,1    3,1
> 1,3    3,3
...
> Those results explain why my image is translating to right and bottom
> when reducing the destination size
>
> I would like to know if it is because of decimation process (nearest
> neighbor) ?
> In this case the returned result should be like this :
> For destination dimension = 2 :
> 0,0    2,0
> 0,2    2,2

Than Hvu Le,

I disagree with your conclusion on how things ought to work.  I think
either decimation system would be consistent with the downsampling.
Currently the sampling takes the center of the destination pixel and
transforms it back to the source image.  In this case it falls on a
boundary between two pixels and which is selected is arbitrary.

If you wish to change the logic the place to look would be the IRasterIO
implementation in the gdal/gcore/rasterio.cpp file - circa line 408
which looks like:

/* -------------------------------------------------------------------- */
/*      Read case                                                       */
/*      Loop over buffer computing source locations.                    */
/* -------------------------------------------------------------------- */
        for( iBufYOff = 0; iBufYOff < nBufYSize; iBufYOff++ )
        {
            size_t   iBufOffset, iSrcOffset;

            dfSrcY = (iBufYOff+0.5) * dfSrcYInc + nYOff;
            iSrcY = (int) dfSrcY;

There is similar logic for the "X" case just below it.  Changing
it here should affect most drivers, though a few handle the
decimation themselves.

Good luck,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent