5
\$\begingroup\$

I am implementing map and reduce - style functions for processing geospatial raster datasets. I would like the map and reduce functions to accept a user-defined function as an input which will be applied to the raster dataset. Currently I am using function pointers - is this a good starting point? Could there be problems scaling this to a large code base, to accept complicated algorithms etc?

Project Aims:

As well as providing the mentioned map and reduce functions for users to pass their own custom functions to, I would eventually like to provide a library of functions which can be passed to map/reduce. As such a library operates on arrays, it may be useful outside of the geospatial domain, so having it decoupled from the map/reduce library is quite useful - users could ignore that and just pass arrays to this 'algorithms' library. I think my initial approach to this would be to wrap selected functions from OpenCV (or another, similar library) which are useful for geospatial raster analysis, before developing other custom functions.

Given the above, I suppose the overall aim could be thought of as a marriage of GDAL and OpenCV while remaining incredibly flexible and extensible.

The Code

So far I have only implemented the map function and the template for a RasterProcess class.

header:

// Include processing functions

/**
 * \brief Definition of a raster processing function.
 *
 * A GALGRasterProcessFn accepts an array of data as input, applies custom logic and writes the output to padfOutArray.
 * Such a function can be passed to GALGRunRasterProcess to apply custom processing to a GDALDataset in chunks and create
 * a new GDALDataset.
 *
 * @param padfInArray The input array of data.
 *
 * @param padfOutArray The output array of data. On first call (via GDALRunRasterProcess) this will be an empty, initialised array,
 *    which should be populated with the result of calculations on padfInArray. In subsequent calls it will contain the result of the
 *    previous window.
 *
 * @param nWindowXSize the actual x size (width) of the read window.
 *
 * @param nWindowYSize the actual y size (height) of the read window. The length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
 *
 * @param pData Process-specific data. This data is passed straight through to the GDALRasterProcessFn and may contain e.g user defined parameters.
 *     The GDALRasterProcessFn definition would define the structure/type of such data.
 *
 * @param pdfNoDataValue The no data value of the dataset
 */

typedef GALGError (*rasterProcessFn)(double *padfInArray, double *padfOutArray,
        int nWindowXSize, int nWindowYSize, void *pData,
        double *pdfInNoDataValue, double *pdfOutNoDataValue);


class GALG_EXPORT RasterProcess {

public:
    RasterProcess();


    /**
     * \brief Apply a raster processing function to each sub-window of a raster.
     *
     * The input raster dataset is read in chunks of nWindowXSize * nWindowYSize and each chunk is passed to the processing
     * function. The output array from the function is written to the destination dataset.
     * An optional 'pixel buffer' can be specified to allow overlaps between successive windows. This is useful for
     * some algorithms, e.g. blob extraction, watershed/stream flow analysis, convolution etc.
     * Process specific data can be passed (e.g. configuration parameters). This data is simply passed straight through to the processing
     * function on each call.
     *
     * @param processFn A GALGRasterProcessFn to apply to each sub window of the raster.
     *
     * @param inputPathStr Path to the source raster dataset from which pixel values are read
     *
     * @param outputPathStr Path to the desired output GeoTiff dataset
     *
     * @param dataObject Process-specific data. This is passed straight through to the GDALRasterProcessFn on each call.
     *
     * @param windowXSize The desired width of each read window. If NULL it defaults to the 'natural' block size of the raster
     *
     * @param windowYSize The desired height of each read window. If NULL it defaults to the 'natural' block size.
     *
     * @param nPixelBuffer A pixel buffer to apply to the read window. The read window is expanded by pnPixelBuffer pixels in all directions such that
     *    each window overlaps by pnPixelBuffer pixels.
     *
     * @param skipHoles If true, will skip processing blocks which contain only no data values and create a sparse geotiff. Only available for geotiff inputs
     *
     * @return a GALGError struct indicating whether the process succeeded.
     */
    GALGError map(rasterProcessFn processFn, const char *inputPathStr,
                const char *outputPathStr, void *dataObject, int *windowXSize,
                int *windowYSize, int *nPixelBuffer, bool skipHoles);

    /**
     * \brief Apply multiple raster processing functions to each sub-window of a raster
     *
     * For each window, the functions defined by the paProcessFn array are called in turn, with the array output of the previous function forming the input
     * to the next function. This allows processing 'toolchains' to be built without having to create intermediate datasets, which can be less efficient in time and space.
     *
     *
     * @param processFnArray An array of GDALRasterProcessFn to apply to each sub window of the raster
     *
     * @param nProcesses The size of paProcessFn
     *
     * @param inputPathStr The path to the source raster dataset from which pixel values are read
     *
     * @param outputPathStr The path to the destination raster dataset to which pixel values are written. Must support RasterIO in write mode.
     *
     * @param dataObjectArray an array of process-specific data objects of size nProcesses. Each data object will be passed to the corresponding GDALRasterProcessFn
     *
     * @param windowXSize The desired width of each read window. If NULL it defaults to the 'natural' block size of the raster
     *
     * @param windowYSize The desired height of each read window. If NULL it defaults to the 'natural' block size.
     *
     * @param nPixelBuffer A pixel buffer to apply to the read window. The read window is expanded by pnPixelBuffer pixels in all directions such that
     *    each window overlaps by pnPixelBuffer pixels.
     *
     *    @param skipHoles If true, will skip processing blocks which contain only no data values and create a sparse geotiff. Only available for geotiff inputs
     *
     * @return a GALGError struct indicating whether the process succeeded.
     */
    GALGError mapMany(rasterProcessFn **processFnArray, int nProcesses,
                    const char *inputPathStr, const char *outputPathStr, void **dataObjectArray,
                    int *windowXSize, int *windowYSize, int *nPixelBuffer, bool skipHoles);


    /**
     * \brief Apply a raster processing 'reduction' function to each sub-window of multiple raster datasets.
     *
     * TODO: Complete
     */
    GALGError reduce(rasterProcessFn processFn, const char **inputPathStrArray,
            const char *outputPathStr, void *dataObject, int *windowXSize,
            int *windowYSize, int *nPixelBuffer, bool skipHoles);


};

implementation:

GALGError createOutputDataset(GDALDataset *srcDataset, const char *outputPathStr, GDALDataset *dstDataset, bool skipHoles) {
    GALGError errResult = { 0, NULL };

    const char *formatStr = "GTiff";
    GDALDriver *gdalDriver;
    gdalDriver = GetGDALDriverManager()->GetDriverByName(formatStr);

    RETURNIF(gdalDriver == NULL, 1, "Could not initialise Geotiff driver")

    char **optionStrArray;
    optionStrArray = CSLSetNameValue(optionStrArray, "TILED", "YES");
    optionStrArray = CSLSetNameValue(optionStrArray, "COMPRESS", "LZW");
    if (skipHoles) {
        optionStrArray = CSLSetNameValue(optionStrArray, "SPARSE_OK", "TRUE");
    }

    dstDataset = gdalDriver->Create(outputPathStr, srcDataset->GetRasterXSize(), srcDataset->GetRasterYSize(),
        srcDataset->GetRasterCount(), srcDataset->GetRasterBand(1)->GetRasterDataType(),
        optionStrArray);

    RETURNIF(dstDataset == NULL, 1, "Could not create output dataset");

    double geotransform[6];
    srcDataset->GetGeoTransform(geotransform);
    dstDataset->SetGeoTransform(geotransform);
    dstDataset->SetProjection(srcDataset->GetProjectionRef());

    GDALRasterBand *srcBand, *dstBand;
    for (int ixBand = 0; ixBand < srcDataset->GetRasterCount(); ++ixBand) {
        srcBand = srcDataset->GetRasterBand(ixBand + 1);
        dstBand = dstDataset->GetRasterBand(ixBand + 1);
        dstBand->SetNoDataValue(srcBand->GetNoDataValue());
    }
    return errResult;
}

RasterProcess::RasterProcess(){

}

GALGError RasterProcess::map(rasterProcessFn processFn, const char *inputPathStr,
        const char *outputPathStr, void *dataObject, int *windowXSize,
        int *windowYSize, int *nPixelBuffer, bool skipHoles) {

    GALGError result = { 0, NULL };
    GDALDataset *srcDataset = NULL, *dstDataset = NULL;

    // Open the input dataset and verify
    srcDataset = (GDALDataset *)GDALOpenEx(inputPathStr, NULL, NULL, NULL, NULL);
    RETURNIF(srcDataset == NULL, 1, "Could not open source dataset");

    // Create output dataset and verify
    result = createOutputDataset(srcDataset, outputPathStr, dstDataset, skipHoles);
    RETURNIF(result.errnum != 0, result.errnum, result.msg);

    // Setup the iterator. If pixelBuffer was passed, we created a buffered iterator,
    // otherwise use a standard BlockIterator
    BlockIterator *iterator = NULL;
    if (nPixelBuffer != NULL) {
        iterator = new BufferedIterator(dstDataset, *nPixelBuffer);
    } else {
        iterator = new BlockIterator(dstDataset);
    }
    RETURNIF(iterator == NULL, 1, "Unable to allocate memory for BlockIterator");
    iterator->setBlockSize(*windowXSize, *windowYSize);

    // Prepare the data buffers
    double *bufInputData = NULL, *bufOutputData = NULL; 
    bufInputData = (double *) VSIMalloc2((size_t) * windowXSize,
            (size_t) * windowYSize);
    bufOutputData = (double *) VSIMalloc2((size_t) * windowXSize,
            (size_t) * windowYSize);
    RETURNIF(bufInputData == NULL || bufOutputData == NULL, 1, "Unable to allocate data arrays");

    int nBands = srcDataset->GetRasterCount();
    int xOff, yOff, xSize, ySize;
    GDALRasterBand *srcBand, *dstBand;
    double inNoDataValue, outNoDataValue;
    int bSuccess;

    // Apply the process function to each sub window of each band
    // in the dataset
    for (int iBand = 0; iBand < nBands; ++iBand) {

        srcBand = srcDataset->GetRasterBand(iBand);
        dstBand = dstDataset->GetRasterBand(iBand);
        inNoDataValue = srcBand->GetNoDataValue(&bSuccess);
        outNoDataValue = dstBand->GetNoDataValue(&bSuccess);

        while (iterator->next(&xSize, &ySize, &xOff, &yOff)) {

            // Read input data. TODO:: Verify function ran
            srcBand->RasterIO(GF_Read, xOff, yOff, xSize, ySize,
                    bufInputData, xSize, ySize, GDT_Float64, 0, 0);

            // Call the process function. TODO: Verify output
            processFn(bufInputData, bufOutputData, xSize, ySize, dataObject,
                    &inNoDataValue, &outNoDataValue);

            // Write out the result. TODO: Verify function ran
            dstBand->RasterIO(GF_Write, xOff, yOff, xSize, ySize,
                    bufOutputData, xSize, ySize, GDT_Float64, 0, 0);
        }
    }

    return result;
}

GALGError RasterProcess::mapMany(rasterProcessFn **processFnArray, int nProcesses,
        const char *inputPathStr, const char *outputPathStr, void **dataObjectArray,
        int *windowXSize, int *windowYSize, int *nPixelBuffer, bool skipHoles) {

    GALGError result = { 1, "Not Implemented" };

    // TODO

    return result;
}

GALGError RasterProcess::reduce(rasterProcessFn processFn, const char **inputPathStrArray,
            const char *outputPathStr, void *dataObject, int *windowXSize,
            int *windowYSize, int *nPixelBuffer, bool skipHoles){
    GALGError result = { 1, "Not Implemented" };

    // TODO

    return result;
}

Obviously, there is still much to implement, but the map function is...functional, and this is the guts of the whole thing.

\$\endgroup\$
2
+50
\$\begingroup\$

To be completely honest, this code snippet appears to be C rather than C++ with a few exceptions.

There is not enough code here to compile it, the included header files are missing. It is not clear due to the missing header files, but it appears that RETURNIF() is a macro. In C++ it would be better to use an inline function than a macro. Macros do not supply any type checking, inline functions supply type checking.

Use References Rather Than Pointers
There are still times when it is proper to use pointers in C++, but generally it is discouraged. In this function declaration

GALGError createOutputDataset(GDALDataset *srcDataset, const char *outputPathStr, GDALDataset *dstDataset, bool skipHoles) {

it appears that pointers are being used to prevent data copying, this same result can be implemented using references rather than pointers. Using references makes the code easier to read and maintain. The same function declaration using references would be

GALGError createOutputDataset(GDALDataset& srcDataset, const char *outputPathStr, GDALDataset& dstDataset, bool skipHoles) {

An example of how the code becomes easier to read and maintain, the following line in createOutputDataset():

    dstDataset = gdalDriver->Create(outputPathStr, srcDataset->GetRasterXSize(), srcDataset->GetRasterYSize(),
        srcDataset->GetRasterCount(), srcDataset->GetRasterBand(1)->GetRasterDataType(),
        optionStrArray);

would be written as :

    dstDataset = gdalDriver.Create(outputPathStr, srcDataset->GetRasterXSize(), srcDataset.GetRasterYSize(),
        srcDataset.GetRasterCount(), srcDataset.GetRasterBand(1).GetRasterDataType(), optionStrArray);

Using references remove the need to pass the address of the object into the function, the call to createOutputDataset() would change from
GALGErrorvalue = createOutputDataset(&srcDataset, *outputPathStr, &dstDataset, skipHoles); to GALGErrorvalue = createOutputDataset(srcDataset, outputPathStr, dstDataset, skipHoles);

This StackOverflow.com question points out the differences between pointers and references and explains when to use them.

Mixing the use of malloc/free with new/delete
The code is mixing new and malloc in the same function

    // Setup the iterator. If pixelBuffer was passed, we created a buffered iterator,
    // otherwise use a standard BlockIterator
    BlockIterator *iterator = NULL;
    if (nPixelBuffer != NULL) {
        iterator = new BufferedIterator(dstDataset, *nPixelBuffer);
    } else {
        iterator = new BlockIterator(dstDataset);
    }
    RETURNIF(iterator == NULL, 1, "Unable to allocate memory for BlockIterator");
    iterator->setBlockSize(*windowXSize, *windowYSize);

    // Prepare the data buffers
    double *bufInputData = NULL, *bufOutputData = NULL; 
    bufInputData = (double *) VSIMalloc2((size_t) * windowXSize,
            (size_t) * windowYSize);
    bufOutputData = (double *) VSIMalloc2((size_t) * windowXSize,
            (size_t) * windowYSize);
    RETURNIF(bufInputData == NULL || bufOutputData == NULL, 1, "Unable to allocate data arrays");

I've seen answers here on Code Review that say that this can cause memory allocation problems, although I can't find this solidly documented on the internet. The new()/delete are prefered because they use the constructors and destructors, and are type safe. Using malloc()/free() is discouraged the constructors and destructors are not used and the allocated memory needs to be cast to the proper type. C type casts are generally discouraged in C++, static_cast and dynamic_cast are prefered.

This StackOverflow.com question goes into detail about why not to use malloc/free in C++, but the basic reason is that new() is type safe, and malloc is not type safe.

This article is also an interesting discussion on malloc()/free() and new()/delete().

|improve this answer|||||
\$\endgroup\$
2
\$\begingroup\$

In modern C++ algorithm is conventionally a template function and the callback is an argument (and automatically deduced template parameter).

// C++03
template <class F>
void traverse(F& cb) { ... }

// C++14
void traverse(auto& cb) { ... }

Such signature lets you pass anything callable, being that a raw pointer, an instance of a class with overloaded () operator, an instance of std::function<...> with compatible signature.

One reason why you might want to accept a more generic callback than just a function pointer is that it is not uncommon for callbacks to have state, which with raw function pointer becomes difficult to implement. Another reason is you might want to be C++1x-friendly and let people use generic lambdas. For instance a prototype code snippet to research your algorithm might be something very localized and easy to type as

map([&](auto&&... args) { std::cout << std::tie(args...) << std::endl; }, ...);

With modern compilers you can pre-instantiate most commonly used templates (say raw pointer case) in the place of algorithm definition to keep object files size to the minimum.

|improve this answer|||||
\$\endgroup\$
  • \$\begingroup\$ I'm curious what does this have to do with the posters code specifically? Can you point out where you might use it? \$\endgroup\$ – pacmaninbw Aug 1 '16 at 11:34
  • \$\begingroup\$ map([](auto&&... args){ std::cout << std::tie(args...) << std::endl; }, ...) to research the algorithm, for instance \$\endgroup\$ – bobah Aug 1 '16 at 13:45
  • 1
    \$\begingroup\$ This help page may help you get some points for your answer if you update your answer codereview.stackexchange.com/help/how-to-answer \$\endgroup\$ – pacmaninbw Aug 1 '16 at 13:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.