// 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);


};

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;
}