[gdal-dev] Proposed Block iterator API

[Even Rouault]

Hi,

I was wondering if there wasn't a need for a block iterator API for code that 
wants to process raster efficiently and iterates over blocks from top to bottom, 
and left to write. Currently you have to create 2 for loops and deal with 
partial blocks at edges (see 
https://trac.osgeo.org/gdal/browser/trunk/gdal/gcore/rasterio.cpp#L3576). 
Nothing really fundamentaly complicated, but I guess people have to often 
write that boring code and I'm sure there a lot of occurences of it inside 
GDAL itself. The API just returns the coordinates of the window, but doesn't 
do the pixel queries (people could use it to operate on a single band, or 
multiple bands)

An example on how to use the API is :

int main(int argc, char* argv[])
{
    GDALAllRegister();
    GDALDatasetH hDS = GDALOpen(argv[1], GA_ReadOnly);
    GDALBlockIteratorH hIter = GDALCreateBlockIterator(
                                                  GDALGetRasterBand(hDS, 1));
    int nXOff, nYOff, nXSize, nYSize;
    while( GDALBlockIteratorGetNext(hIter, &nXOff, &nYOff, &nXSize, &nYSize) )
    {
        double dfProgressStart, dfProgressEnd;
        dfProgressStart = GDALBlockIteratorGetProgress(hIter, &dfProgressEnd);
        printf("nXOff=%d nYOff=%d nXSize=%d nYSize=%d [%.2f->%.2f]\n",
               nXOff, nYOff, nXSize, nYSize, dfProgressStart, dfProgressEnd);
    }
    GDALDestroyBlockIterator(hIter);
    GDALClose(hDS);
    return 0;
}

And if you run it on a 20x20 raster with blocks of size 16x16, it will 
display:

nXOff=0 nYOff=0 nXSize=16 nYSize=16 [0.00->64.00]
nXOff=16 nYOff=0 nXSize=4 nYSize=16 [64.00->80.00]
nXOff=0 nYOff=16 nXSize=16 nYSize=4 [80.00->96.00]
nXOff=16 nYOff=16 nXSize=4 nYSize=4 [96.00->100.00]:

Thoughts ?

Best regards,

Even



PS: below the API and its implementation

#include <gdal_priv.h>

/* API */
CPL_C_START
typedef struct GDALBlockIterator *GDALBlockIteratorH;

GDALBlockIteratorH CPL_DLL GDALCreateBlockIterator(GDALRasterBandH hBand);
void CPL_DLL GDALDestroyBlockIterator(GDALBlockIteratorH hIter);
int CPL_DLL GDALBlockIteratorGetMaxXSize(GDALBlockIteratorH hIter);
int CPL_DLL GDALBlockIteratorGetMaxYSize(GDALBlockIteratorH hIter);
void GDALBlockIteratorSetBlockSize(GDALBlockIterator* poIter, int nBlockWidth, 
int nBlockHeight);
int CPL_DLL GDALBlockIteratorGetNext(GDALBlockIteratorH hIter,
                                     int* pnXOff,
                                     int* pnYOff,
                                     int* pnXSize,
                                     int* pnYSize);
double CPL_DLL GDALBlockIteratorGetProgress(GDALBlockIteratorH hIter, double* 
pdfNextProgress);
CPL_C_END

/* Impl */
struct GDALBlockIterator
{
    int nXSize;
    int nYSize;
    int nXOff;
    int nYOff;
    int nBlockWidth;
    int nBlockHeight;
    bool bFinished;
};
typedef struct GDALBlockIterator GDALBlockIterator;


GDALBlockIterator* GDALCreateBlockIterator(GDALRasterBandH hBand)
{
    GDALRasterBand* poBand = (GDALRasterBand*)hBand;
    GDALBlockIterator* poIter = new GDALBlockIterator;
    poIter->nXSize = poBand->GetXSize();
    poIter->nYSize = poBand->GetYSize();
    poBand->GetBlockSize( &(poIter->nBlockWidth), &(poIter->nBlockHeight) );
    poIter->nXOff = -1;
    poIter->nYOff = 0;
    poIter->bFinished = false;
    return poIter;
}

void GDALDestroyBlockIterator(GDALBlockIteratorH hIter)
{
    delete hIter;
}

int GDALBlockIteratorGetMaxXSize(GDALBlockIterator* poIter)
{
    return poIter->nBlockWidth;
}

int GDALBlockIteratorGetMaxYSize(GDALBlockIterator* poIter)
{
    return poIter->nBlockHeight;
}

void GDALBlockIteratorSetBlockSize(GDALBlockIterator* poIter, int nBlockWidth, 
int nBlockHeight)
{
    poIter->nBlockWidth = nBlockWidth;
    poIter->nBlockHeight = nBlockHeight;
}

int GDALBlockIteratorGetNext(GDALBlockIterator* poIter,
                             int* pnXOff,
                             int* pnYOff,
                             int* pnXSize,
                             int* pnYSize)
{
    if( poIter->bFinished )
        return FALSE;

    if( poIter->nXOff < 0 )
        poIter->nXOff = 0;
    else
        poIter->nXOff += poIter->nBlockWidth;
    if( poIter->nXOff >= poIter->nXSize )
    {
        poIter->nXOff = 0;
        poIter->nYOff += poIter->nBlockHeight;
        if( poIter->nYOff >= poIter->nYSize )
        {
            poIter->bFinished = true;
            return FALSE;
        }
    }

    if( pnXOff )
        *pnXOff = poIter->nXOff;
    if( pnYOff )
        *pnYOff = poIter->nYOff;
    if( pnXSize )
    {
        *pnXSize = poIter->nBlockWidth;
        if( poIter->nXOff + poIter->nBlockWidth >= poIter->nXSize )
        {
            *pnXSize = poIter->nXSize - poIter->nXOff;
        }
    }
    if( pnYSize )
    {
        *pnYSize = poIter->nBlockHeight;
        if(poIter->nYOff + poIter->nBlockHeight >= poIter->nYSize )
        {
            *pnYSize = poIter->nYSize - poIter->nYOff;
        }
    }
    return TRUE;
}

double GDALBlockIteratorGetProgress(GDALBlockIterator* poIter, double* 
pdfNextProgress)
{
    if( poIter->bFinished )
    {
        if( pdfNextProgress != NULL )
            *pdfNextProgress = 100.0;
        return 100.0;
    }

    int nBlockHeight = poIter->nBlockHeight;
    if( poIter->nYOff + nBlockHeight > poIter->nYSize )
        nBlockHeight = poIter->nYSize - poIter->nYOff;
    double dfProcessedPixels = (double)poIter->nYOff * poIter->nXSize + poIter-
>nXOff * nBlockHeight;
    double dfRet = 100.0 * dfProcessedPixels / poIter->nXSize / poIter-
>nYSize;
    if( pdfNextProgress != NULL )
    {
        GDALBlockIterator sIterDup = *poIter;
        GDALBlockIteratorGetNext(&sIterDup, NULL, NULL, NULL, NULL);
        *pdfNextProgress = GDALBlockIteratorGetProgress(&sIterDup, NULL);
    }
    return dfRet;
}

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com