/******************************************************************************
 *
 * Project:  High Performance Image Reprojector
 * Purpose:  Implementation of high level convenience APIs for warper.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "gdalwarper.h"

#include <stdlib.h>
#include <string.h>

#include <algorithm>
#include <cmath>
#include <limits>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_float.h"
#include "cpl_mask.h"
#include "cpl_minixml.h"
#include "cpl_progress.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_priv.h"
#include "ogr_api.h"
#include "ogr_core.h"
#include "vrtdataset.h"  // for VRTSerializeNoData

#if (defined(__x86_64) || defined(_M_X64))
#include <emmintrin.h>
#endif

/************************************************************************/
/*                         GDALReprojectImage()                         */
/************************************************************************/

/**
 * Reproject image.
 *
 * This is a convenience function utilizing the GDALWarpOperation class to
 * reproject an image from a source to a destination.  In particular, this
 * function takes care of establishing the transformation function to
 * implement the reprojection, and will default a variety of other
 * warp options.
 *
 * Nodata values set on destination dataset are taken into account.
 *
 * No metadata, projection info, or color tables are transferred
 * to the output file. Source overviews are not considered.
 *
 * For more advanced warping capabilities, consider using GDALWarp().
 *
 * @param hSrcDS the source image file.
 * @param pszSrcWKT the source projection.  If NULL the source projection
 * is read from from hSrcDS.
 * @param hDstDS the destination image file.
 * @param pszDstWKT the destination projection.  If NULL the destination
 * projection will be read from hDstDS.
 * @param eResampleAlg the type of resampling to use.
 * @param dfWarpMemoryLimit the amount of memory (in bytes) that the warp
 * API is allowed to use for caching.  This is in addition to the memory
 * already allocated to the GDAL caching (as per GDALSetCacheMax()).  May be
 * 0.0 to use default memory settings.
 * @param dfMaxError maximum error measured in input pixels that is allowed
 * in approximating the transformation (0.0 for exact calculations).
 * @param pfnProgress a GDALProgressFunc() compatible callback function for
 * reporting progress or NULL.
 * @param pProgressArg argument to be passed to pfnProgress.  May be NULL.
 * @param psOptions warp options, normally NULL.
 *
 * @return CE_None on success or CE_Failure if something goes wrong.
 * @see GDALWarp()
 */

CPLErr CPL_STDCALL GDALReprojectImage(
    GDALDatasetH hSrcDS, const char *pszSrcWKT, GDALDatasetH hDstDS,
    const char *pszDstWKT, GDALResampleAlg eResampleAlg,
    CPL_UNUSED double dfWarpMemoryLimit, double dfMaxError,
    GDALProgressFunc pfnProgress, void *pProgressArg,
    GDALWarpOptions *psOptions)

{
    /* -------------------------------------------------------------------- */
    /*      Setup a reprojection based transformer.                         */
    /* -------------------------------------------------------------------- */
    void *hTransformArg = GDALCreateGenImgProjTransformer(
        hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, TRUE, 1000.0, 0);

    if (hTransformArg == nullptr)
        return CE_Failure;

    /* -------------------------------------------------------------------- */
    /*      Create a copy of the user provided options, or a defaulted      */
    /*      options structure.                                              */
    /* -------------------------------------------------------------------- */
    GDALWarpOptions *psWOptions = psOptions == nullptr
                                      ? GDALCreateWarpOptions()
                                      : GDALCloneWarpOptions(psOptions);

    psWOptions->eResampleAlg = eResampleAlg;

    /* -------------------------------------------------------------------- */
    /*      Set transform.                                                  */
    /* -------------------------------------------------------------------- */
    if (dfMaxError > 0.0)
    {
        psWOptions->pTransformerArg = GDALCreateApproxTransformer(
            GDALGenImgProjTransform, hTransformArg, dfMaxError);

        psWOptions->pfnTransformer = GDALApproxTransform;
    }
    else
    {
        psWOptions->pfnTransformer = GDALGenImgProjTransform;
        psWOptions->pTransformerArg = hTransformArg;
    }

    /* -------------------------------------------------------------------- */
    /*      Set file and band mapping.                                      */
    /* -------------------------------------------------------------------- */
    psWOptions->hSrcDS = hSrcDS;
    psWOptions->hDstDS = hDstDS;

    int nSrcBands = GDALGetRasterCount(hSrcDS);
    {
        GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, nSrcBands);
        if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
        {
            psWOptions->nSrcAlphaBand = nSrcBands;
            nSrcBands--;
        }
    }

    int nDstBands = GDALGetRasterCount(hDstDS);
    {
        GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, nDstBands);
        if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
        {
            psWOptions->nDstAlphaBand = nDstBands;
            nDstBands--;
        }
    }

    GDALWarpInitDefaultBandMapping(psWOptions, std::min(nSrcBands, nDstBands));

    /* -------------------------------------------------------------------- */
    /*      Set source nodata values if the source dataset seems to have    */
    /*      any. Same for target nodata values                              */
    /* -------------------------------------------------------------------- */
    for (int iBand = 0; iBand < psWOptions->nBandCount; iBand++)
    {
        GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, iBand + 1);

        int bGotNoData = FALSE;
        double dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData);
        if (bGotNoData)
        {
            GDALWarpInitSrcNoDataReal(psWOptions, -1.1e20);
            psWOptions->padfSrcNoDataReal[iBand] = dfNoDataValue;
        }

        // Deal with target band.
        hBand = GDALGetRasterBand(hDstDS, iBand + 1);

        dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData);
        if (bGotNoData)
        {
            GDALWarpInitDstNoDataReal(psWOptions, -1.1e20);
            psWOptions->padfDstNoDataReal[iBand] = dfNoDataValue;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Set the progress function.                                      */
    /* -------------------------------------------------------------------- */
    if (pfnProgress != nullptr)
    {
        psWOptions->pfnProgress = pfnProgress;
        psWOptions->pProgressArg = pProgressArg;
    }

    /* -------------------------------------------------------------------- */
    /*      Create a warp options based on the options.                     */
    /* -------------------------------------------------------------------- */
    GDALWarpOperation oWarper;
    CPLErr eErr = oWarper.Initialize(psWOptions);

    if (eErr == CE_None)
        eErr = oWarper.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS),
                                         GDALGetRasterYSize(hDstDS));

    /* -------------------------------------------------------------------- */
    /*      Cleanup.                                                        */
    /* -------------------------------------------------------------------- */
    GDALDestroyGenImgProjTransformer(hTransformArg);

    if (dfMaxError > 0.0)
        GDALDestroyApproxTransformer(psWOptions->pTransformerArg);

    GDALDestroyWarpOptions(psWOptions);

    return eErr;
}

/************************************************************************/
/*                    GDALCreateAndReprojectImage()                     */
/*                                                                      */
/*      This is a "quicky" reprojection API.                            */
/************************************************************************/

/** Reproject an image and create the target reprojected image */
CPLErr CPL_STDCALL GDALCreateAndReprojectImage(
    GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstFilename,
    const char *pszDstWKT, GDALDriverH hDstDriver, char **papszCreateOptions,
    GDALResampleAlg eResampleAlg, double dfWarpMemoryLimit, double dfMaxError,
    GDALProgressFunc pfnProgress, void *pProgressArg,
    GDALWarpOptions *psOptions)

{
    VALIDATE_POINTER1(hSrcDS, "GDALCreateAndReprojectImage", CE_Failure);

    /* -------------------------------------------------------------------- */
    /*      Default a few parameters.                                       */
    /* -------------------------------------------------------------------- */
    if (hDstDriver == nullptr)
    {
        hDstDriver = GDALGetDriverByName("GTiff");
        if (hDstDriver == nullptr)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "GDALCreateAndReprojectImage needs GTiff driver");
            return CE_Failure;
        }
    }

    if (pszSrcWKT == nullptr)
        pszSrcWKT = GDALGetProjectionRef(hSrcDS);

    if (pszDstWKT == nullptr)
        pszDstWKT = pszSrcWKT;

    /* -------------------------------------------------------------------- */
    /*      Create a transformation object from the source to               */
    /*      destination coordinate system.                                  */
    /* -------------------------------------------------------------------- */
    void *hTransformArg = GDALCreateGenImgProjTransformer(
        hSrcDS, pszSrcWKT, nullptr, pszDstWKT, TRUE, 1000.0, 0);

    if (hTransformArg == nullptr)
        return CE_Failure;

    /* -------------------------------------------------------------------- */
    /*      Get approximate output definition.                              */
    /* -------------------------------------------------------------------- */
    double adfDstGeoTransform[6] = {};
    int nPixels = 0;
    int nLines = 0;

    if (GDALSuggestedWarpOutput(hSrcDS, GDALGenImgProjTransform, hTransformArg,
                                adfDstGeoTransform, &nPixels,
                                &nLines) != CE_None)
        return CE_Failure;

    GDALDestroyGenImgProjTransformer(hTransformArg);

    /* -------------------------------------------------------------------- */
    /*      Create the output file.                                         */
    /* -------------------------------------------------------------------- */
    GDALDatasetH hDstDS = GDALCreate(
        hDstDriver, pszDstFilename, nPixels, nLines, GDALGetRasterCount(hSrcDS),
        GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1)),
        papszCreateOptions);

    if (hDstDS == nullptr)
        return CE_Failure;

    /* -------------------------------------------------------------------- */
    /*      Write out the projection definition.                            */
    /* -------------------------------------------------------------------- */
    GDALSetProjection(hDstDS, pszDstWKT);
    GDALSetGeoTransform(hDstDS, adfDstGeoTransform);

    /* -------------------------------------------------------------------- */
    /*      Perform the reprojection.                                       */
    /* -------------------------------------------------------------------- */
    CPLErr eErr = GDALReprojectImage(
        hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, eResampleAlg, dfWarpMemoryLimit,
        dfMaxError, pfnProgress, pProgressArg, psOptions);

    GDALClose(hDstDS);

    return eErr;
}

/************************************************************************/
/*                       GDALWarpNoDataMaskerT()                        */
/************************************************************************/

template <class T>
static CPLErr GDALWarpNoDataMaskerT(const double *padfNoData, size_t nPixels,
                                    const T *pData, GUInt32 *panValidityMask,
                                    int *pbOutAllValid)
{
    // Nothing to do if value is out of range.
    if (padfNoData[0] < cpl::NumericLimits<T>::min() ||
        padfNoData[0] > cpl::NumericLimits<T>::max() + 0.000001 ||
        padfNoData[1] != 0.0)
    {
        *pbOutAllValid = TRUE;
        return CE_None;
    }

    const int nNoData = static_cast<int>(floor(padfNoData[0] + 0.000001));
    int bAllValid = TRUE;
    for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
    {
        if (pData[iOffset] == nNoData)
        {
            bAllValid = FALSE;
            CPLMaskClear(panValidityMask, iOffset);
        }
    }
    *pbOutAllValid = bAllValid;

    return CE_None;
}

/************************************************************************/
/*                        GDALWarpNoDataMasker()                        */
/*                                                                      */
/*      GDALMaskFunc for establishing a validity mask for a source      */
/*      band based on a provided NODATA value.                          */
/************************************************************************/

CPLErr GDALWarpNoDataMasker(void *pMaskFuncArg, int nBandCount,
                            GDALDataType eType, int /* nXOff */,
                            int /* nYOff */, int nXSize, int nYSize,
                            GByte **ppImageData, int bMaskIsFloat,
                            void *pValidityMask, int *pbOutAllValid)

{
    const double *padfNoData = static_cast<double *>(pMaskFuncArg);
    GUInt32 *panValidityMask = static_cast<GUInt32 *>(pValidityMask);
    const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;

    *pbOutAllValid = FALSE;

    if (nBandCount != 1 || bMaskIsFloat)
    {
        CPLError(
            CE_Failure, CPLE_AppDefined,
            "Invalid nBandCount or bMaskIsFloat argument in SourceNoDataMask");
        return CE_Failure;
    }

    CPLErr eErr = CE_None;

    switch (eType)
    {
        case GDT_Byte:
            return GDALWarpNoDataMaskerT(padfNoData, nPixels,
                                         *ppImageData,  // Already a GByte *.
                                         panValidityMask, pbOutAllValid);

        case GDT_Int16:
            return GDALWarpNoDataMaskerT(
                padfNoData, nPixels, reinterpret_cast<GInt16 *>(*ppImageData),
                panValidityMask, pbOutAllValid);

        case GDT_UInt16:
            return GDALWarpNoDataMaskerT(
                padfNoData, nPixels, reinterpret_cast<GUInt16 *>(*ppImageData),
                panValidityMask, pbOutAllValid);

        case GDT_Float32:
        {
            const float fNoData = static_cast<float>(padfNoData[0]);
            const float *pafData = reinterpret_cast<float *>(*ppImageData);
            const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(fNoData));

            // Nothing to do if value is out of range.
            if (padfNoData[1] != 0.0)
            {
                *pbOutAllValid = TRUE;
                return CE_None;
            }

            int bAllValid = TRUE;
            for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
            {
                float fVal = pafData[iOffset];
                if ((bIsNoDataNan && std::isnan(fVal)) ||
                    (!bIsNoDataNan && ARE_REAL_EQUAL(fVal, fNoData)))
                {
                    bAllValid = FALSE;
                    CPLMaskClear(panValidityMask, iOffset);
                }
            }
            *pbOutAllValid = bAllValid;
        }
        break;

        case GDT_Float64:
        {
            const double dfNoData = padfNoData[0];
            const double *padfData = reinterpret_cast<double *>(*ppImageData);
            const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(dfNoData));

            // Nothing to do if value is out of range.
            if (padfNoData[1] != 0.0)
            {
                *pbOutAllValid = TRUE;
                return CE_None;
            }

            int bAllValid = TRUE;
            for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
            {
                double dfVal = padfData[iOffset];
                if ((bIsNoDataNan && std::isnan(dfVal)) ||
                    (!bIsNoDataNan && ARE_REAL_EQUAL(dfVal, dfNoData)))
                {
                    bAllValid = FALSE;
                    CPLMaskClear(panValidityMask, iOffset);
                }
            }
            *pbOutAllValid = bAllValid;
        }
        break;

        default:
        {
            const int nWordSize = GDALGetDataTypeSizeBytes(eType);

            const bool bIsNoDataRealNan =
                CPL_TO_BOOL(std::isnan(padfNoData[0]));

            eErr = CE_Failure;
            double *padfWrk = static_cast<double *>(
                VSI_MALLOC2_VERBOSE(nXSize, sizeof(double) * 2));
            if (padfWrk)
            {
                eErr = CE_None;
                bool bAllValid = true;
                for (int iLine = 0; iLine < nYSize; iLine++)
                {
                    GDALCopyWords((*ppImageData) + nWordSize * iLine * nXSize,
                                  eType, nWordSize, padfWrk, GDT_CFloat64, 16,
                                  nXSize);

                    for (int iPixel = 0; iPixel < nXSize; ++iPixel)
                    {
                        if (((bIsNoDataRealNan &&
                              std::isnan(padfWrk[iPixel * 2])) ||
                             (!bIsNoDataRealNan &&
                              ARE_REAL_EQUAL(padfWrk[iPixel * 2],
                                             padfNoData[0]))))
                        {
                            size_t iOffset =
                                iPixel + static_cast<size_t>(iLine) * nXSize;

                            bAllValid = false;
                            CPLMaskClear(panValidityMask, iOffset);
                        }
                    }
                }
                *pbOutAllValid = bAllValid;

                VSIFree(padfWrk);
            }
        }
        break;
    }

    return eErr;
}

/************************************************************************/
/*                       GDALWarpSrcAlphaMasker()                       */
/*                                                                      */
/*      GDALMaskFunc for reading source simple 8bit alpha mask          */
/*      information and building a floating point density mask from     */
/*      it.                                                             */
/************************************************************************/

CPLErr GDALWarpSrcAlphaMasker(void *pMaskFuncArg, int /* nBandCount */,
                              GDALDataType /* eType */, int nXOff, int nYOff,
                              int nXSize, int nYSize, GByte ** /*ppImageData */,
                              int bMaskIsFloat, void *pValidityMask,
                              int *pbOutAllOpaque)

{
    GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
    float *pafMask = static_cast<float *>(pValidityMask);
    *pbOutAllOpaque = FALSE;
    const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;

    /* -------------------------------------------------------------------- */
    /*      Do some minimal checking.                                       */
    /* -------------------------------------------------------------------- */
    if (!bMaskIsFloat)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    if (psWO == nullptr || psWO->nSrcAlphaBand < 1)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    /* -------------------------------------------------------------------- */
    /*      Read the alpha band.                                            */
    /* -------------------------------------------------------------------- */
    GDALRasterBandH hAlphaBand =
        GDALGetRasterBand(psWO->hSrcDS, psWO->nSrcAlphaBand);
    if (hAlphaBand == nullptr)
        return CE_Failure;

    // Rescale.
    const float inv_alpha_max = static_cast<float>(
        1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions,
                                           "SRC_ALPHA_MAX", "255")));
    bool bOutAllOpaque = true;

    size_t iPixel = 0;
    CPLErr eErr;

#if (defined(__x86_64) || defined(_M_X64))
    GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
    // Make sure that pafMask is at least 8-byte aligned, which should
    // normally be always the case if being a ptr returned by malloc().
    if ((eDT == GDT_Byte || eDT == GDT_UInt16) && CPL_IS_ALIGNED(pafMask, 8))
    {
        // Read data.
        eErr = GDALRasterIOEx(
            hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask, nXSize,
            nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
            static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);

        if (eErr != CE_None)
            return eErr;

        // Make sure we have the correct alignment before doing SSE
        // On Linux x86_64, the alignment should be always correct due
        // the alignment of malloc() being 16 byte.
        const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff;
        if (!CPL_IS_ALIGNED(pafMask, 16))
        {
            pafMask[iPixel] =
                (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
                inv_alpha_max;
            if (pafMask[iPixel] >= 1.0f)
                pafMask[iPixel] = 1.0f;
            else
                bOutAllOpaque = false;
            iPixel++;
        }
        CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
        const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max);
        const float one_single = 1.0f;
        const __m128 xmm_one = _mm_load1_ps(&one_single);
        const __m128i xmm_i_mask = _mm_set1_epi32(mask);
        __m128 xmmMaskNonOpaque0 = _mm_setzero_ps();
        __m128 xmmMaskNonOpaque1 = _mm_setzero_ps();
        __m128 xmmMaskNonOpaque2 = _mm_setzero_ps();
        for (; iPixel + 6 * 4 - 1 < nPixels; iPixel += 6 * 4)
        {
            __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 0))));
            __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 1))));
            __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 2))));
            __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 3))));
            __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 4))));
            __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128(
                xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                pafMask + iPixel + 4 * 5))));
            xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max);
            xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max);
            xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max);
            xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max);
            xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max);
            xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max);
            xmmMaskNonOpaque0 =
                _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask0, xmm_one));
            xmmMaskNonOpaque1 =
                _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask1, xmm_one));
            xmmMaskNonOpaque2 =
                _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask2, xmm_one));
            xmmMaskNonOpaque0 =
                _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask3, xmm_one));
            xmmMaskNonOpaque1 =
                _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask4, xmm_one));
            xmmMaskNonOpaque2 =
                _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask5, xmm_one));
            xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one);
            xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one);
            xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one);
            xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one);
            xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one);
            xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one);
            _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0);
            _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1);
            _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2);
            _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3);
            _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4);
            _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5);
        }
        if (_mm_movemask_ps(
                _mm_or_ps(_mm_or_ps(xmmMaskNonOpaque0, xmmMaskNonOpaque1),
                          xmmMaskNonOpaque2)))
        {
            bOutAllOpaque = false;
        }
        for (; iPixel < nPixels; iPixel++)
        {
            pafMask[iPixel] =
                (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
                inv_alpha_max;
            if (pafMask[iPixel] >= 1.0f)
                pafMask[iPixel] = 1.0f;
            else
                bOutAllOpaque = false;
        }
    }
    else
#endif
    {
        // Read data.
        eErr = GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
                            pafMask, nXSize, nYSize, GDT_Float32, 0, 0);

        if (eErr != CE_None)
            return eErr;

        // TODO(rouault): Is loop unrolling by hand (r34564) actually helpful?
        for (; iPixel + 3 < nPixels; iPixel += 4)
        {
            pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
            if (pafMask[iPixel] >= 1.0f)
                pafMask[iPixel] = 1.0f;
            else
                bOutAllOpaque = false;
            pafMask[iPixel + 1] = pafMask[iPixel + 1] * inv_alpha_max;
            if (pafMask[iPixel + 1] >= 1.0f)
                pafMask[iPixel + 1] = 1.0f;
            else
                bOutAllOpaque = false;
            pafMask[iPixel + 2] = pafMask[iPixel + 2] * inv_alpha_max;
            if (pafMask[iPixel + 2] >= 1.0f)
                pafMask[iPixel + 2] = 1.0f;
            else
                bOutAllOpaque = false;
            pafMask[iPixel + 3] = pafMask[iPixel + 3] * inv_alpha_max;
            if (pafMask[iPixel + 3] >= 1.0f)
                pafMask[iPixel + 3] = 1.0f;
            else
                bOutAllOpaque = false;
        }

        for (; iPixel < nPixels; iPixel++)
        {
            pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
            if (pafMask[iPixel] >= 1.0f)
                pafMask[iPixel] = 1.0f;
            else
                bOutAllOpaque = false;
        }
    }

    *pbOutAllOpaque = bOutAllOpaque;

    return CE_None;
}

/************************************************************************/
/*                       GDALWarpSrcMaskMasker()                        */
/*                                                                      */
/*      GDALMaskFunc for reading source simple 8bit validity mask       */
/*      information and building a one bit validity mask.               */
/************************************************************************/

CPLErr GDALWarpSrcMaskMasker(void *pMaskFuncArg, int /* nBandCount */,
                             GDALDataType /* eType */, int nXOff, int nYOff,
                             int nXSize, int nYSize, GByte ** /*ppImageData */,
                             int bMaskIsFloat, void *pValidityMask)

{
    GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
    GUInt32 *panMask = static_cast<GUInt32 *>(pValidityMask);

    /* -------------------------------------------------------------------- */
    /*      Do some minimal checking.                                       */
    /* -------------------------------------------------------------------- */
    if (bMaskIsFloat)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    if (psWO == nullptr)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    /* -------------------------------------------------------------------- */
    /*      Allocate a temporary buffer to read mask byte data into.        */
    /* -------------------------------------------------------------------- */
    GByte *pabySrcMask =
        static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nYSize));
    if (pabySrcMask == nullptr)
    {
        return CE_Failure;
    }

    /* -------------------------------------------------------------------- */
    /*      Fetch our mask band.                                            */
    /* -------------------------------------------------------------------- */
    GDALRasterBandH hMaskBand = nullptr;
    GDALRasterBandH hSrcBand =
        GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[0]);
    if (hSrcBand != nullptr)
        hMaskBand = GDALGetMaskBand(hSrcBand);

    if (hMaskBand == nullptr)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    /* -------------------------------------------------------------------- */
    /*      Read the mask band.                                             */
    /* -------------------------------------------------------------------- */
    CPLErr eErr = GDALRasterIO(hMaskBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
                               pabySrcMask, nXSize, nYSize, GDT_Byte, 0, 0);

    if (eErr != CE_None)
    {
        CPLFree(pabySrcMask);
        return eErr;
    }

    /* -------------------------------------------------------------------- */
    /*      Pack into 1 bit per pixel for validity.                         */
    /* -------------------------------------------------------------------- */
    const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;
    for (size_t iPixel = 0; iPixel < nPixels; iPixel++)
    {
        if (pabySrcMask[iPixel] == 0)
            CPLMaskClear(panMask, iPixel);
    }

    CPLFree(pabySrcMask);

    return CE_None;
}

/************************************************************************/
/*                       GDALWarpDstAlphaMasker()                       */
/*                                                                      */
/*      GDALMaskFunc for reading or writing the destination simple      */
/*      8bit alpha mask information and building a floating point       */
/*      density mask from it.   Note, writing is distinguished          */
/*      negative bandcount.                                             */
/************************************************************************/

CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
                              CPL_UNUSED GDALDataType /* eType */, int nXOff,
                              int nYOff, int nXSize, int nYSize,
                              GByte ** /*ppImageData */, int bMaskIsFloat,
                              void *pValidityMask)
{
    /* -------------------------------------------------------------------- */
    /*      Do some minimal checking.                                       */
    /* -------------------------------------------------------------------- */
    if (!bMaskIsFloat)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
    if (psWO == nullptr || psWO->nDstAlphaBand < 1)
    {
        CPLAssert(false);
        return CE_Failure;
    }

    float *pafMask = static_cast<float *>(pValidityMask);
    const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;

    GDALRasterBandH hAlphaBand =
        GDALGetRasterBand(psWO->hDstDS, psWO->nDstAlphaBand);
    if (hAlphaBand == nullptr)
        return CE_Failure;

    size_t iPixel = 0;

    /* -------------------------------------------------------------------- */
    /*      Read alpha case.                                                */
    /* -------------------------------------------------------------------- */
    if (nBandCount >= 0)
    {
        const char *pszInitDest =
            CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST");

        // Special logic for destinations being initialized on-the-fly.
        if (pszInitDest != nullptr)
        {
            memset(pafMask, 0, nPixels * sizeof(float));
            return CE_None;
        }

        // Rescale.
        const float inv_alpha_max = static_cast<float>(
            1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions,
                                               "DST_ALPHA_MAX", "255")));

#if (defined(__x86_64) || defined(_M_X64))
        const GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
        // Make sure that pafMask is at least 8-byte aligned, which should
        // normally be always the case if being a ptr returned by malloc().
        if ((eDT == GDT_Byte || eDT == GDT_UInt16) &&
            CPL_IS_ALIGNED(pafMask, 8))
        {
            // Read data.
            const CPLErr eErr = GDALRasterIOEx(
                hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask,
                nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
                static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);

            if (eErr != CE_None)
                return eErr;

            // Make sure we have the correct alignment before doing SSE
            // On Linux x86_64, the alignment should be always correct due
            // the alignment of malloc() being 16 byte.
            const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff;
            if (!CPL_IS_ALIGNED(pafMask, 16))
            {
                pafMask[iPixel] =
                    (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
                    inv_alpha_max;
                pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
                iPixel++;
            }
            CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
            const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max);
            const float one_single = 1.0f;
            const __m128 xmm_one = _mm_load1_ps(&one_single);
            const __m128i xmm_i_mask = _mm_set1_epi32(mask);
            for (; iPixel + 31 < nPixels; iPixel += 32)
            {
                __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 0))));
                __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 1))));
                __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 2))));
                __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 3))));
                __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 4))));
                __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 5))));
                __m128 xmm_mask6 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 6))));
                __m128 xmm_mask7 = _mm_cvtepi32_ps(_mm_and_si128(
                    xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
                                    pafMask + iPixel + 4 * 7))));
                xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max);
                xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max);
                xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max);
                xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max);
                xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max);
                xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max);
                xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_inverse_alpha_max);
                xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_inverse_alpha_max);
                xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one);
                xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one);
                xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one);
                xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one);
                xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one);
                xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one);
                xmm_mask6 = _mm_min_ps(xmm_mask6, xmm_one);
                xmm_mask7 = _mm_min_ps(xmm_mask7, xmm_one);
                _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0);
                _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1);
                _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2);
                _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3);
                _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4);
                _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5);
                _mm_store_ps(pafMask + iPixel + 4 * 6, xmm_mask6);
                _mm_store_ps(pafMask + iPixel + 4 * 7, xmm_mask7);
            }
            for (; iPixel < nPixels; iPixel++)
            {
                pafMask[iPixel] =
                    (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
                    inv_alpha_max;
                pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
            }
        }
        else
#endif
        {
            // Read data.
            const CPLErr eErr =
                GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
                             pafMask, nXSize, nYSize, GDT_Float32, 0, 0);

            if (eErr != CE_None)
                return eErr;

            for (; iPixel < nPixels; iPixel++)
            {
                pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
                pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
            }
        }

        return CE_None;
    }

    /* -------------------------------------------------------------------- */
    /*      Write alpha case.                                               */
    /* -------------------------------------------------------------------- */
    else
    {
        GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
        const float cst_alpha_max =
            static_cast<float>(CPLAtof(CSLFetchNameValueDef(
                psWO->papszWarpOptions, "DST_ALPHA_MAX", "255"))) +
            ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16 ||
              eDT == GDT_Int32 || eDT == GDT_UInt32)
                 ? 0.1f
                 : 0.0f);

        CPLErr eErr = CE_None;

#if (defined(__x86_64) || defined(_M_X64))
        // Make sure that pafMask is at least 8-byte aligned, which should
        // normally be always the case if being a ptr returned by malloc()
        if ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16) &&
            CPL_IS_ALIGNED(pafMask, 8))
        {
            // Make sure we have the correct alignment before doing SSE
            // On Linux x86_64, the alignment should be always correct due
            // the alignment of malloc() being 16 byte
            if (!CPL_IS_ALIGNED(pafMask, 16))
            {
                reinterpret_cast<int *>(pafMask)[iPixel] =
                    static_cast<int>(pafMask[iPixel] * cst_alpha_max);
                iPixel++;
            }
            CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
            const __m128 xmm_alpha_max = _mm_load1_ps(&cst_alpha_max);
            for (; iPixel + 31 < nPixels; iPixel += 32)
            {
                __m128 xmm_mask0 = _mm_load_ps(pafMask + iPixel + 4 * 0);
                __m128 xmm_mask1 = _mm_load_ps(pafMask + iPixel + 4 * 1);
                __m128 xmm_mask2 = _mm_load_ps(pafMask + iPixel + 4 * 2);
                __m128 xmm_mask3 = _mm_load_ps(pafMask + iPixel + 4 * 3);
                __m128 xmm_mask4 = _mm_load_ps(pafMask + iPixel + 4 * 4);
                __m128 xmm_mask5 = _mm_load_ps(pafMask + iPixel + 4 * 5);
                __m128 xmm_mask6 = _mm_load_ps(pafMask + iPixel + 4 * 6);
                __m128 xmm_mask7 = _mm_load_ps(pafMask + iPixel + 4 * 7);
                xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_alpha_max);
                xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_alpha_max);
                xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_alpha_max);
                xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_alpha_max);
                xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_alpha_max);
                xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_alpha_max);
                xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_alpha_max);
                xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_alpha_max);
                // Truncate to int.
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 0),
                    _mm_cvttps_epi32(xmm_mask0));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 1),
                    _mm_cvttps_epi32(xmm_mask1));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 2),
                    _mm_cvttps_epi32(xmm_mask2));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 3),
                    _mm_cvttps_epi32(xmm_mask3));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 4),
                    _mm_cvttps_epi32(xmm_mask4));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 5),
                    _mm_cvttps_epi32(xmm_mask5));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 6),
                    _mm_cvttps_epi32(xmm_mask6));
                _mm_store_si128(
                    reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 7),
                    _mm_cvttps_epi32(xmm_mask7));
            }
            for (; iPixel < nPixels; iPixel++)
                reinterpret_cast<int *>(pafMask)[iPixel] =
                    static_cast<int>(pafMask[iPixel] * cst_alpha_max);

            // Write data.
            // Assumes little endianness here.
            eErr = GDALRasterIOEx(
                hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize, pafMask,
                nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
                static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);
        }
        else
#endif
        {
            for (; iPixel + 3 < nPixels; iPixel += 4)
            {
                pafMask[iPixel + 0] = static_cast<float>(
                    static_cast<int>(pafMask[iPixel + 0] * cst_alpha_max));
                pafMask[iPixel + 1] = static_cast<float>(
                    static_cast<int>(pafMask[iPixel + 1] * cst_alpha_max));
                pafMask[iPixel + 2] = static_cast<float>(
                    static_cast<int>(pafMask[iPixel + 2] * cst_alpha_max));
                pafMask[iPixel + 3] = static_cast<float>(
                    static_cast<int>(pafMask[iPixel + 3] * cst_alpha_max));
            }
            for (; iPixel < nPixels; iPixel++)
                pafMask[iPixel] = static_cast<float>(
                    static_cast<int>(pafMask[iPixel] * cst_alpha_max));

            // Write data.

            eErr =
                GDALRasterIO(hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize,
                             pafMask, nXSize, nYSize, GDT_Float32, 0, 0);
        }
        return eErr;
    }
}

/************************************************************************/
/*                      GDALWarpGetOptionList()                         */
/************************************************************************/

/** Return a XML string describing options accepted by
 * GDALWarpOptions::papszWarpOptions.
 *
 * @since 3.11
 */
const char *GDALWarpGetOptionList(void)
{
    return "<OptionList>"
           "<Option name='INIT_DEST' type='string' description='"
           "Numeric value or NO_DATA. This option forces the destination image "
           "to be initialized to the indicated value (for all bands) "
           "or indicates that it should be initialized to the NO_DATA value in "
           "padfDstNoDataReal/padfDstNoDataImag. If this value is not set the "
           "destination image will be read and overlaid.'/>"
           "<Option name='WRITE_FLUSH' type='boolean' description='"
           "This option forces a flush to disk of data after "
           "each chunk is processed. In some cases this helps ensure a serial "
           " writing of the output data otherwise a block of data may be "
           "written to disk each time a block of data is read for the input "
           "buffer resulting in a lot of extra seeking around the disk, and "
           "reduced IO throughput.' default='NO'/>"
           "<Option name='SKIP_NOSOURCE' type='boolean' description='"
           "Skip all processing for chunks for which there is no corresponding "
           "input data. This will disable initializing the destination "
           "(INIT_DEST) and all other processing, and so should be used "
           "carefully.  Mostly useful to short circuit a lot of extra work "
           "in mosaicing situations. gdalwarp will automatically enable this "
           "option when it is assumed to be safe to do so.' default='NO'/>"
#ifdef undocumented
           "<Option name='ERROR_OUT_IF_EMPTY_SOURCE_WINDOW' type='boolean' "
           "description='By default, if the source window corresponding to the "
           "current target window fails to be determined due to reprojection "
           "errors, the warping fails. Setting this option to NO prevent such "
           "failure from happening. The warped VRT mechanism automatically "
           "sets it to NO.'/>"
#endif
           "<Option name='UNIFIED_SRC_NODATA' type='string-select' "
           "description='"
           "This setting determines how to take into account nodata values "
           "when there are several input bands. Consult "
           "GDALWarpOptions::papszWarpOptions documentation for more details.'>"
           "  <Value>AUTO</Value>"
           "  <Value>PARTIAL</Value>"
           "  <Value>YES</Value>"
           "  <Value>NO</Value>"
           "</Option>"
           "<Option name='CUTLINE' type='string' description='"
           "This may contain the WKT geometry for a cutline.  It will be "
           "converted into a geometry by GDALWarpOperation::Initialize() and "
           "assigned to the GDALWarpOptions hCutline field. The coordinates "
           "must be expressed in source pixel/line coordinates. Note: this is "
           "different from the assumptions made for the -cutline option "
           "of the gdalwarp utility !'/>"
           "<Option name='CUTLINE_BLEND_DIST' type='float' description='"
           "This may be set with a distance in pixels which will be assigned "
           "to the dfCutlineBlendDist field in the GDALWarpOptions.'/>"
           "<Option name='CUTLINE_ALL_TOUCHED' type='boolean' description='"
           "This may be set to TRUE to enable ALL_TOUCHED mode when "
           "rasterizing cutline polygons. This is useful to ensure that that "
           "all pixels overlapping the cutline polygon will be selected, not "
           "just those whose center point falls within the polygon.' "
           "default='NO'/>"
           "<Option name='XSCALE' type='float' description='"
           "Ratio expressing the resampling factor (number of destination "
           "pixels per source pixel) along the target horizontal axis. The "
           "scale is used to determine the number of source pixels along the "
           "x-axis that are considered by the resampling algorithm. "
           "Equals to one for no resampling, below one for downsampling "
           "and above one for upsampling. This is automatically computed, "
           "for each processing chunk, and may thus vary among them, depending "
           "on the shape of output regions vs input regions. Such variations "
           "can be undesired in some situations. If the resampling factor "
           "can be considered as constant over the warped area, setting a "
           "constant value can lead to more reproducible pixel output.'/>"
           "<Option name='YSCALE' type='float' description='"
           "Same as XSCALE, but along the horizontal axis.'/>"
           "<Option name='OPTIMIZE_SIZE' type='boolean' description='"
           "This defaults to FALSE, but may be set to TRUE typically when "
           "writing to a compressed dataset (GeoTIFF with COMPRESS creation "
           "option set for example) for achieving a smaller file size. This "
           "is achieved by writing at once data aligned on full blocks of the "
           "target dataset, which avoids partial writes of compressed blocks "
           "and lost space when they are rewritten at the end of the file. "
           "However sticking to target block size may cause major processing "
           "slowdown for some particular reprojections. OPTIMIZE_SIZE mode "
           "is automatically enabled when it is safe to do so. "
           "As this parameter influences the shape of warping chunk, and by "
           "default the XSCALE and YSCALE parameters are computed per warping "
           "chunk, this parameter may influence the pixel output.' "
           "default='NO'/>"
           "<Option name='NUM_THREADS' type='string' description='"
           "Can be set to a numeric value or ALL_CPUS to set the number of "
           "threads to use to parallelize the computation part of the warping. "
           "If not set, computation will be done in a single thread..'/>"
           "<Option name='STREAMABLE_OUTPUT' type='boolean' description='"
           "This defaults to FALSE, but may be set to TRUE typically when "
           "writing to a streamed file. The gdalwarp utility automatically "
           "sets this option when writing to /vsistdout/ or a named pipe "
           "(on Unix). This option has performance impacts for some "
           "reprojections. Note: band interleaved output is "
           "not currently supported by the warping algorithm in a streamable "
           "compatible way.' default='NO'/>"
           "<Option name='SRC_COORD_PRECISION' type='float' description='"
           "Advanced setting. This defaults to 0, to indicate that no rounding "
           "of computing source image coordinates corresponding to the target "
           "image must be done. If greater than 0 (and typically below 1), "
           "this value, expressed in pixel, will be used to round computed "
           "source image coordinates. The purpose of this option is to make "
           "the results of warping with the approximated transformer more "
           "reproducible and not sensitive to changes in warping memory size. "
           "To achieve that, SRC_COORD_PRECISION must be at least 10 times "
           "greater than the error threshold. The higher the "
           "SRC_COORD_PRECISION/error_threshold ratio, the higher the "
           "performance will be, since exact reprojections must statistically "
           "be done with a frequency of "
           "4*error_threshold/SRC_COORD_PRECISION.' default='0'/>"
           "<Option name='SRC_ALPHA_MAX' type='float' description='"
           "Maximum value for the alpha band of the source dataset. If the "
           "value is not set and the alpha band has a NBITS metadata item, "
           "it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if the "
           "value is not set and the alpha band is of type UInt16 "
           "(resp Int16), 65535 (resp 32767) is used. "
           "Otherwise, 255 is used.'/>"
           "<Option name='DST_ALPHA_MAX' type='float' description='"
           "Maximum value for the alpha band of the destination dataset. "
           "If the value is not set and the alpha band has a NBITS metadata "
           "item, it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if "
           "the value is not set and the alpha band is of type UInt16 "
           "(resp Int16), 65535 (resp 32767) is used. "
           "Otherwise, 255 is used.'/>"
           "<Option name='SAMPLE_GRID' type='boolean' description='"
           "Setting this option to YES will force the sampling to "
           "include internal points as well as edge points which can be "
           "important if the transformation is esoteric inside out, or if "
           "large sections of the destination image are not transformable into "
           "the source coordinate system.' default='NO'/>"
           "<Option name='SAMPLE_STEPS' type='string' description='"
           "Modifies the density of the sampling grid. Increasing this can "
           "increase the computational cost, but improves the accuracy with "
           "which the source region is computed. This can be set to ALL to "
           "mean to sample along all edge points of the destination region "
           "(if SAMPLE_GRID=NO or not specified), or all points of the "
           "destination region if SAMPLE_GRID=YES.' default='21'/>"
           "<Option name='SOURCE_EXTRA' type='int' description='"
           "This is a number of extra pixels added around the source "
           "window for a given request, and by default it is 1 to take care "
           "of rounding error. Setting this larger will increase the amount of "
           "data that needs to be read, but can avoid missing source data.' "
           "default='1'/>"
           "<Option name='APPLY_VERTICAL_SHIFT' type='boolean' description='"
           "Force the use of vertical shift. This option is generally not "
           "necessary, except when using an explicit coordinate transformation "
           "(COORDINATE_OPERATION), and not specifying an explicit source and "
           "target SRS.'/>"
           "<Option name='MULT_FACTOR_VERTICAL_SHIFT' type='float' "
           "description='"
           "Multiplication factor for the vertical shift' default='1.0'/>"
           "<Option name='EXCLUDED_VALUES' type='string' "
           "description='"
           "Comma-separated tuple of values (thus typically \"R,G,B\"), that "
           "are ignored as contributing source pixels during resampling. "
           "The number of values in the tuple must be the same as the number "
           "of bands, excluding the alpha band. Several tuples of excluded "
           "values may be specified using the \"(R1,G1,B2),(R2,G2,B2)\" syntax."
           " Only taken into account by Average currently. This concept is a "
           "bit similar to nodata/alpha, but the main difference is that "
           "pixels matching one of the excluded value tuples are still "
           "considered as valid, when determining the target pixel "
           "validity/density.'/>"
           "<Option name='EXCLUDED_VALUES_PCT_THRESHOLD' type='float' "
           "min='0' max='100' description='"
           "Minimum percentage of source pixels that must be set at one of "
           "the EXCLUDED_VALUES to cause the excluded value, that is in "
           "majority among source pixels, to be used as the target pixel "
           "value. Only taken into account by Average currently.' "
           "default='50'/>"
           "<Option name='NODATA_VALUES_PCT_THRESHOLD' type='float' "
           "min='0' max='100' description='"
           "Minimum percentage of source pixels that must be at nodata (or "
           "alpha=0 or any other way to express transparent pixel) to cause "
           "the target pixel value to not be set. Default value is 100 (%), "
           "which means that a target pixel is not set only if all "
           "contributing source pixels are not set. Note that "
           "NODATA_VALUES_PCT_THRESHOLD is taken into account before "
           "EXCLUDED_VALUES_PCT_THRESHOLD. Only taken into account by Average "
           "currently.' default='100'/>"
           "<Option name='MODE_TIES' type='string-select' "
           "description='"
           "Strategy to use when breaking ties with MODE resampling. "
           "By default, the first value encountered will be used. "
           "Alternatively, the minimum or maximum value can be selected.' "
           "default='FIRST'>"
           "  <Value>FIRST</Value>"
           "  <Value>MIN</Value>"
           "  <Value>MAX</Value>"
           "</Option>"
           "</OptionList>";
}

/************************************************************************/
/* ==================================================================== */
/*                           GDALWarpOptions                            */
/* ==================================================================== */
/************************************************************************/

/**
 * \var char **GDALWarpOptions::papszWarpOptions;
 *
 * A string list of additional options controlling the warp operation in
 * name=value format.  A suitable string list can be prepared with
 * CSLSetNameValue().
 *
 * The available options can also be retrieved programmatically with
 * GDALWarpGetOptionList().
 *
 * The following values are currently supported:
 * <ul>
 * <li>INIT_DEST=[value] or INIT_DEST=NO_DATA: This option forces the
 * destination image to be initialized to the indicated value (for all bands)
 * or indicates that it should be initialized to the NO_DATA value in
 * padfDstNoDataReal/padfDstNoDataImag. If this value isn't set the
 * destination image will be read and overlaid.</li>
 *
 * <li>WRITE_FLUSH=YES/NO: This option forces a flush to disk of data after
 * each chunk is processed. In some cases this helps ensure a serial
 * writing of the output data otherwise a block of data may be written to disk
 * each time a block of data is read for the input buffer resulting in a lot
 * of extra seeking around the disk, and reduced IO throughput. The default
 * is NO.</li>
 *
 * <li>SKIP_NOSOURCE=YES/NO: Skip all processing for chunks for which there
 * is no corresponding input data. This will disable initializing the
 * destination (INIT_DEST) and all other processing, and so should be used
 * carefully.  Mostly useful to short circuit a lot of extra work in mosaicing
 * situations. gdalwarp will automatically enable this
 * option when it is assumed to be safe to do so.</li>
 *
 * <li>UNIFIED_SRC_NODATA=YES/NO/PARTIAL: This setting determines
 * how to take into account nodata values when there are several input bands.
 * <ul>
 * <li>When YES, all bands are considered as nodata if and only if, all bands
 *     match the corresponding nodata values.
 *     Note: UNIFIED_SRC_NODATA=YES is set by default, when called from gdalwarp
 * / GDALWarp() with an explicit -srcnodata setting.
 *
 *     Example with nodata values at (1, 2, 3) and target alpha band requested.
 *     <ul>
 *     <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li>
 *     <li>input pixel = (1, 2, 127) ==> output pixel = (1, 2, 127, 255)</li>
 *     </ul>
 * </li>
 * <li>When NO, nodata masking values is considered independently for each band.
 *     A potential target alpha band will always be valid if there are multiple
 *     bands.
 *
 *     Example with nodata values at (1, 2, 3) and target alpha band requested.
 *     <ul>
 *     <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 255)</li>
 *     <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li>
 *     </ul>
 *
 *     Note: NO was the default behavior before GDAL 3.3.2
 * </li>
 * <li>When PARTIAL, or not specified at all (default behavior),
 *     nodata masking values is considered independently for each band.
 *     But, and this is the difference with NO, if for a given pixel, it
 *     evaluates to the nodata value of each band, the target pixel is
 *     considered as globally invalid, which impacts the value of a potential
 *     target alpha band.
 *
 *     Note: PARTIAL is new to GDAL 3.3.2 and should not be used with
 *     earlier versions. The default behavior of GDAL < 3.3.2 was NO.
 *
 *     Example with nodata values at (1, 2, 3) and target alpha band requested.
 *     <ul>
 *     <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li>
 *     <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li>
 *     </ul>
 * </li>
 * </ul>
 * </li>
 *
 * <li>CUTLINE: This may contain the WKT geometry for a cutline.  It will
 * be converted into a geometry by GDALWarpOperation::Initialize() and assigned
 * to the GDALWarpOptions hCutline field. The coordinates must be expressed
 * in source pixel/line coordinates. Note: this is different from the
 * assumptions made for the -cutline option of the gdalwarp utility !</li>
 *
 * <li>CUTLINE_BLEND_DIST: This may be set with a distance in pixels which
 * will be assigned to the dfCutlineBlendDist field in the GDALWarpOptions.</li>
 *
 * <li>CUTLINE_ALL_TOUCHED: This defaults to FALSE, but may be set to TRUE
 * to enable ALL_TOUCHED mode when rasterizing cutline polygons.  This is
 * useful to ensure that that all pixels overlapping the cutline polygon
 * will be selected, not just those whose center point falls within the
 * polygon.</li>
 *
 * <li>XSCALE: Ratio expressing the resampling factor (number of destination
 * pixels per source pixel) along the target horizontal axis.
 * The scale is used to determine the number of source pixels along the x-axis
 * that are considered by the resampling algorithm.
 * Equals to one for no resampling, below one for downsampling
 * and above one for upsampling. This is automatically computed, for each
 * processing chunk, and may thus vary among them, depending on the
 * shape of output regions vs input regions. Such variations can be undesired
 * in some situations. If the resampling factor can be considered as constant
 * over the warped area, setting a constant value can lead to more reproducible
 * pixel output.</li>
 *
 * <li>YSCALE: Same as XSCALE, but along the horizontal axis.</li>
 *
 * <li>OPTIMIZE_SIZE: This defaults to FALSE, but may be set to TRUE
 * typically when writing to a compressed dataset (GeoTIFF with
 * COMPRESS creation option set for example) for achieving a smaller
 * file size. This is achieved by writing at once data aligned on full
 * blocks of the target dataset, which avoids partial writes of
 * compressed blocks and lost space when they are rewritten at the end
 * of the file. However sticking to target block size may cause major
 * processing slowdown for some particular reprojections. Starting
 * with GDAL 3.8, OPTIMIZE_SIZE mode is automatically enabled when it is safe
 * to do so.
 * As this parameter influences the shape of warping chunk, and by default the
 * XSCALE and YSCALE parameters are computed per warping chunk, this parameter may
 * influence the pixel output.
 * </li>
 *
 * <li>NUM_THREADS: (GDAL >= 1.10) Can be set to a numeric value or ALL_CPUS to
 * set the number of threads to use to parallelize the computation part of the
 * warping. If not set, computation will be done in a single thread.</li>
 *
 * <li>STREAMABLE_OUTPUT: This defaults to FALSE, but may
 * be set to TRUE typically when writing to a streamed file. The
 * gdalwarp utility automatically sets this option when writing to
 * /vsistdout/ or a named pipe (on Unix).  This option has performance
 * impacts for some reprojections.  Note: band interleaved output is
 * not currently supported by the warping algorithm in a streamable
 * compatible way.</li>
 *
 * <li>SRC_COORD_PRECISION: Advanced setting. This
 * defaults to 0, to indicate that no rounding of computing source
 * image coordinates corresponding to the target image must be
 * done. If greater than 0 (and typically below 1), this value,
 * expressed in pixel, will be used to round computed source image
 * coordinates. The purpose of this option is to make the results of
 * warping with the approximated transformer more reproducible and not
 * sensitive to changes in warping memory size. To achieve that,
 * SRC_COORD_PRECISION must be at least 10 times greater than the
 * error threshold. The higher the SRC_COORD_PRECISION/error_threshold
 * ratio, the higher the performance will be, since exact
 * reprojections must statistically be done with a frequency of
 * 4*error_threshold/SRC_COORD_PRECISION.</li>
 *
 * <li>SRC_ALPHA_MAX: Maximum value for the alpha band of the
 * source dataset. If the value is not set and the alpha band has a NBITS
 * metadata item, it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if the
 * value is not set and the alpha band is of type UInt16 (resp Int16), 65535
 * (resp 32767) is used. Otherwise, 255 is used.</li>
 *
 * <li>DST_ALPHA_MAX: Maximum value for the alpha band of the
 * destination dataset. If the value is not set and the alpha band has a NBITS
 * metadata item, it is used to set DST_ALPHA_MAX = 2^NBITS-1. Otherwise, if the
 * value is not set and the alpha band is of type UInt16 (resp Int16), 65535
 * (resp 32767) is used. Otherwise, 255 is used.</li>
 * </ul>
 *
 * Normally when computing the source raster data to
 * load to generate a particular output area, the warper samples transforms
 * 21 points along each edge of the destination region back onto the source
 * file, and uses this to compute a bounding window on the source image that
 * is sufficient.  Depending on the transformation in effect, the source
 * window may be a bit too small, or even missing large areas.  Problem
 * situations are those where the transformation is very non-linear or
 * "inside out".  Examples are transforming from WGS84 to Polar Stereographic
 * for areas around the pole, or transformations where some of the image is
 * untransformable.  The following options provide some additional control
 * to deal with errors in computing the source window:
 * <ul>
 *
 * <li>SAMPLE_GRID=YES/NO: Setting this option to YES will force the sampling to
 * include internal points as well as edge points which can be important if
 * the transformation is esoteric inside out, or if large sections of the
 * destination image are not transformable into the source coordinate
 * system.</li>
 *
 * <li>SAMPLE_STEPS: Modifies the density of the sampling grid.  The default
 * number of steps is 21.   Increasing this can increase the computational
 * cost, but improves the accuracy with which the source region is
 * computed.
 * Starting with GDAL 3.7, this can be set to ALL to mean to sample
 * along all edge points of the destination region (if SAMPLE_GRID=NO or not
 * specified), or all points of the destination region if SAMPLE_GRID=YES.</li>
 *
 * <li>SOURCE_EXTRA: This is a number of extra pixels added around the source
 * window for a given request, and by default it is 1 to take care of rounding
 * error.  Setting this larger will increase the amount of data that needs to
 * be read, but can avoid missing source data.</li>
 * <li>APPLY_VERTICAL_SHIFT=YES/NO: Force the use of vertical shift.
 * This option is generally not necessary, except when using an explicit
 * coordinate transformation (COORDINATE_OPERATION), and not specifying
 * an explicit source and target SRS.</li>
 * <li>MULT_FACTOR_VERTICAL_SHIFT: Multiplication factor for the vertical
 * shift. Default 1.0</li>
 *
 * <li>EXCLUDED_VALUES: (GDAL >= 3.9) Comma-separated tuple of values
 * (thus typically "R,G,B"), that are ignored as contributing source
 * pixels during resampling. The number of values in the tuple must be the same
 * as the number of bands, excluding the alpha band.
 * Several tuples of excluded values may be specified using the
 * "(R1,G1,B2),(R2,G2,B2)" syntax.
 * Only taken into account by Average currently.
 * This concept is a bit similar to nodata/alpha, but the main difference is
 * that pixels matching one of the excluded value tuples are still considered
 * as valid, when determining the target pixel validity/density.
 * </li>
 *
 * <li>EXCLUDED_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage
 * of source pixels that must be set at one of the EXCLUDED_VALUES to cause
 * the excluded value, that is in majority among source pixels, to be used as the
 * target pixel value. Default value is 50 (%).
 * Only taken into account by Average currently.</li>
 *
 * <li>NODATA_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage
 * of source pixels that must be at nodata (or alpha=0 or any other way to express
 * transparent pixel) to cause the target pixel value to not be set. Default
 * value is 100 (%), which means that a target pixel is not set only if all
 * contributing source pixels are not set.
 * Note that NODATA_VALUES_PCT_THRESHOLD is taken into account before
 * EXCLUDED_VALUES_PCT_THRESHOLD.
 * Only taken into account by Average currently.</li>
 *
 * <li>MODE_TIES=FIRST/MIN/MAX: (GDAL >= 3.11) Strategy to use when breaking
 * ties with MODE resampling. By default, the first value encountered will be used.
 * Alternatively, the minimum or maximum value can be selected.</li>
 *
 * </ul>
 */

/************************************************************************/
/*                       GDALCreateWarpOptions()                        */
/************************************************************************/

/** Create a warp options structure.
 *
 * Must be deallocated with GDALDestroyWarpOptions()
 */
GDALWarpOptions *CPL_STDCALL GDALCreateWarpOptions()

{
    GDALWarpOptions *psOptions =
        static_cast<GDALWarpOptions *>(CPLCalloc(sizeof(GDALWarpOptions), 1));

    psOptions->nBandCount = 0;
    psOptions->eResampleAlg = GRA_NearestNeighbour;
    psOptions->pfnProgress = GDALDummyProgress;
    psOptions->eWorkingDataType = GDT_Unknown;
    psOptions->eTieStrategy = GWKTS_First;

    return psOptions;
}

/************************************************************************/
/*                       GDALDestroyWarpOptions()                       */
/************************************************************************/

/** Destroy a warp options structure. */
void CPL_STDCALL GDALDestroyWarpOptions(GDALWarpOptions *psOptions)

{
    if (psOptions == nullptr)
        return;

    CSLDestroy(psOptions->papszWarpOptions);
    CPLFree(psOptions->panSrcBands);
    CPLFree(psOptions->panDstBands);
    CPLFree(psOptions->padfSrcNoDataReal);
    CPLFree(psOptions->padfSrcNoDataImag);
    CPLFree(psOptions->padfDstNoDataReal);
    CPLFree(psOptions->padfDstNoDataImag);
    CPLFree(psOptions->papfnSrcPerBandValidityMaskFunc);
    CPLFree(psOptions->papSrcPerBandValidityMaskFuncArg);

    if (psOptions->hCutline != nullptr)
        delete static_cast<OGRGeometry *>(psOptions->hCutline);

    CPLFree(psOptions);
}

#define COPY_MEM(target, type, count)                                          \
    do                                                                         \
    {                                                                          \
        if ((psSrcOptions->target) != nullptr && (count) != 0)                 \
        {                                                                      \
            (psDstOptions->target) =                                           \
                static_cast<type *>(CPLMalloc(sizeof(type) * (count)));        \
            memcpy((psDstOptions->target), (psSrcOptions->target),             \
                   sizeof(type) * (count));                                    \
        }                                                                      \
        else                                                                   \
            (psDstOptions->target) = nullptr;                                  \
    } while (false)

/************************************************************************/
/*                        GDALCloneWarpOptions()                        */
/************************************************************************/

/** Clone a warp options structure.
 *
 * Must be deallocated with GDALDestroyWarpOptions()
 */
GDALWarpOptions *CPL_STDCALL
GDALCloneWarpOptions(const GDALWarpOptions *psSrcOptions)

{
    GDALWarpOptions *psDstOptions = GDALCreateWarpOptions();

    memcpy(psDstOptions, psSrcOptions, sizeof(GDALWarpOptions));

    if (psSrcOptions->papszWarpOptions != nullptr)
        psDstOptions->papszWarpOptions =
            CSLDuplicate(psSrcOptions->papszWarpOptions);

    COPY_MEM(panSrcBands, int, psSrcOptions->nBandCount);
    COPY_MEM(panDstBands, int, psSrcOptions->nBandCount);
    COPY_MEM(padfSrcNoDataReal, double, psSrcOptions->nBandCount);
    COPY_MEM(padfSrcNoDataImag, double, psSrcOptions->nBandCount);
    COPY_MEM(padfDstNoDataReal, double, psSrcOptions->nBandCount);
    COPY_MEM(padfDstNoDataImag, double, psSrcOptions->nBandCount);
    // cppcheck-suppress pointerSize
    COPY_MEM(papfnSrcPerBandValidityMaskFunc, GDALMaskFunc,
             psSrcOptions->nBandCount);
    psDstOptions->papSrcPerBandValidityMaskFuncArg = nullptr;

    if (psSrcOptions->hCutline != nullptr)
        psDstOptions->hCutline =
            OGR_G_Clone(static_cast<OGRGeometryH>(psSrcOptions->hCutline));
    psDstOptions->dfCutlineBlendDist = psSrcOptions->dfCutlineBlendDist;

    return psDstOptions;
}

namespace
{
void InitNoData(int nBandCount, double **ppdNoDataReal, double dDataReal)
{
    if (nBandCount <= 0)
    {
        return;
    }
    if (*ppdNoDataReal != nullptr)
    {
        return;
    }

    *ppdNoDataReal =
        static_cast<double *>(CPLMalloc(sizeof(double) * nBandCount));

    for (int i = 0; i < nBandCount; ++i)
    {
        (*ppdNoDataReal)[i] = dDataReal;
    }
}
}  // namespace

/************************************************************************/
/*                      GDALWarpInitDstNoDataReal()                     */
/************************************************************************/

/**
 * \brief Initialize padfDstNoDataReal with specified value.
 *
 * @param psOptionsIn options to initialize.
 * @param dNoDataReal value to initialize to.
 *
 */
void CPL_STDCALL GDALWarpInitDstNoDataReal(GDALWarpOptions *psOptionsIn,
                                           double dNoDataReal)
{
    VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataReal");
    InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataReal,
               dNoDataReal);
}

/************************************************************************/
/*                      GDALWarpInitSrcNoDataReal()                     */
/************************************************************************/

/**
 * \brief Initialize padfSrcNoDataReal with specified value.
 *
 * @param psOptionsIn options to initialize.
 * @param dNoDataReal value to initialize to.
 *
 */
void CPL_STDCALL GDALWarpInitSrcNoDataReal(GDALWarpOptions *psOptionsIn,
                                           double dNoDataReal)
{
    VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataReal");
    InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataReal,
               dNoDataReal);
}

/************************************************************************/
/*                      GDALWarpInitNoDataReal()                        */
/************************************************************************/

/**
 * \brief Initialize padfSrcNoDataReal and padfDstNoDataReal with specified
 * value.
 *
 * @param psOptionsIn options to initialize.
 * @param dNoDataReal value to initialize to.
 *
 */
void CPL_STDCALL GDALWarpInitNoDataReal(GDALWarpOptions *psOptionsIn,
                                        double dNoDataReal)
{
    GDALWarpInitDstNoDataReal(psOptionsIn, dNoDataReal);
    GDALWarpInitSrcNoDataReal(psOptionsIn, dNoDataReal);
}

/************************************************************************/
/*                      GDALWarpInitDstNoDataImag()                     */
/************************************************************************/

/**
 * \brief Initialize padfDstNoDataImag  with specified value.
 *
 * @param psOptionsIn options to initialize.
 * @param dNoDataImag value to initialize to.
 *
 */
void CPL_STDCALL GDALWarpInitDstNoDataImag(GDALWarpOptions *psOptionsIn,
                                           double dNoDataImag)
{
    VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataImag");
    InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataImag,
               dNoDataImag);
}

/************************************************************************/
/*                      GDALWarpInitSrcNoDataImag()                     */
/************************************************************************/

/**
 * \brief Initialize padfSrcNoDataImag  with specified value.
 *
 * @param psOptionsIn options to initialize.
 * @param dNoDataImag value to initialize to.
 *
 */
void CPL_STDCALL GDALWarpInitSrcNoDataImag(GDALWarpOptions *psOptionsIn,
                                           double dNoDataImag)
{
    VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataImag");
    InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataImag,
               dNoDataImag);
}

/************************************************************************/
/*                      GDALWarpResolveWorkingDataType()                */
/************************************************************************/

/**
 * \brief If the working data type is unknown, this method will determine
 *  a valid working data type to support the data in the src and dest
 *  data sets and any noData values.
 *
 * @param psOptions options to initialize.
 *
 */
void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions)
{
    if (psOptions == nullptr)
    {
        return;
    }
    /* -------------------------------------------------------------------- */
    /*      If no working data type was provided, set one now.              */
    /*                                                                      */
    /*      Ensure that the working data type can encapsulate any value     */
    /*      in the target, source, and the no data for either.              */
    /* -------------------------------------------------------------------- */
    if (psOptions->eWorkingDataType != GDT_Unknown)
    {
        return;
    }

    psOptions->eWorkingDataType = GDT_Byte;

    // If none of the provided input nodata values can be represented in the
    // data type of the corresponding source band, ignore them.
    if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal)
    {
        int nCountInvalidSrcNoDataReal = 0;
        for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
        {
            GDALRasterBandH hSrcBand = GDALGetRasterBand(
                psOptions->hSrcDS, psOptions->panSrcBands[iBand]);

            if (hSrcBand &&
                !GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand],
                                    GDALGetRasterDataType(hSrcBand)))
            {
                nCountInvalidSrcNoDataReal++;
            }
        }
        if (nCountInvalidSrcNoDataReal == psOptions->nBandCount)
        {
            CPLFree(psOptions->padfSrcNoDataReal);
            psOptions->padfSrcNoDataReal = nullptr;
            CPLFree(psOptions->padfSrcNoDataImag);
            psOptions->padfSrcNoDataImag = nullptr;
        }
    }

    for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
    {
        if (psOptions->hDstDS != nullptr)
        {
            GDALRasterBandH hDstBand = GDALGetRasterBand(
                psOptions->hDstDS, psOptions->panDstBands[iBand]);

            if (hDstBand != nullptr)
            {
                psOptions->eWorkingDataType =
                    GDALDataTypeUnion(psOptions->eWorkingDataType,
                                      GDALGetRasterDataType(hDstBand));
            }
        }

        if (psOptions->hSrcDS != nullptr)
        {
            GDALRasterBandH hSrcBand = GDALGetRasterBand(
                psOptions->hSrcDS, psOptions->panSrcBands[iBand]);

            if (hSrcBand != nullptr)
            {
                psOptions->eWorkingDataType =
                    GDALDataTypeUnion(psOptions->eWorkingDataType,
                                      GDALGetRasterDataType(hSrcBand));
            }
        }

        if (psOptions->padfSrcNoDataReal != nullptr)
        {
            psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
                psOptions->eWorkingDataType,
                psOptions->padfSrcNoDataReal[iBand], false);
        }

        if (psOptions->padfSrcNoDataImag != nullptr &&
            psOptions->padfSrcNoDataImag[iBand] != 0.0)
        {
            psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
                psOptions->eWorkingDataType,
                psOptions->padfSrcNoDataImag[iBand], true);
        }

        if (psOptions->padfDstNoDataReal != nullptr)
        {
            psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
                psOptions->eWorkingDataType,
                psOptions->padfDstNoDataReal[iBand], false);
        }

        if (psOptions->padfDstNoDataImag != nullptr &&
            psOptions->padfDstNoDataImag[iBand] != 0.0)
        {
            psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
                psOptions->eWorkingDataType,
                psOptions->padfDstNoDataImag[iBand], true);
        }
    }

    const bool bApplyVerticalShift = CPLFetchBool(
        psOptions->papszWarpOptions, "APPLY_VERTICAL_SHIFT", false);
    if (bApplyVerticalShift &&
        GDALDataTypeIsInteger(psOptions->eWorkingDataType))
    {
        const double dfMultFactorVerticalShift = CPLAtof(CSLFetchNameValueDef(
            psOptions->papszWarpOptions, "MULT_FACTOR_VERTICAL_SHIFT", "1.0"));
        if (dfMultFactorVerticalShift != 1)
        {
            psOptions->eWorkingDataType =
                GDALDataTypeUnion(psOptions->eWorkingDataType, GDT_Float32);
        }
    }
}

/************************************************************************/
/*                      GDALWarpInitDefaultBandMapping()                */
/************************************************************************/

/**
 * \brief Init src and dst band mappings such that Bands[i] = i+1
 *  for nBandCount
 *  Does nothing if psOptionsIn->nBandCount is non-zero.
 *
 * @param psOptionsIn options to initialize.
 * @param nBandCount bands to initialize for.
 *
 */
void CPL_STDCALL GDALWarpInitDefaultBandMapping(GDALWarpOptions *psOptionsIn,
                                                int nBandCount)
{
    if (psOptionsIn->nBandCount != 0)
    {
        return;
    }

    psOptionsIn->nBandCount = nBandCount;

    psOptionsIn->panSrcBands =
        static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount));
    psOptionsIn->panDstBands =
        static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount));

    for (int i = 0; i < psOptionsIn->nBandCount; i++)
    {
        psOptionsIn->panSrcBands[i] = i + 1;
        psOptionsIn->panDstBands[i] = i + 1;
    }
}

/************************************************************************/
/*                      GDALSerializeWarpOptions()                      */
/************************************************************************/

CPLXMLNode *CPL_STDCALL GDALSerializeWarpOptions(const GDALWarpOptions *psWO)

{
    /* -------------------------------------------------------------------- */
    /*      Create root.                                                    */
    /* -------------------------------------------------------------------- */
    CPLXMLNode *psTree =
        CPLCreateXMLNode(nullptr, CXT_Element, "GDALWarpOptions");

    /* -------------------------------------------------------------------- */
    /*      WarpMemoryLimit                                                 */
    /* -------------------------------------------------------------------- */
    CPLCreateXMLElementAndValue(
        psTree, "WarpMemoryLimit",
        CPLString().Printf("%g", psWO->dfWarpMemoryLimit));

    /* -------------------------------------------------------------------- */
    /*      ResampleAlg                                                     */
    /* -------------------------------------------------------------------- */
    const char *pszAlgName = nullptr;

    if (psWO->eResampleAlg == GRA_NearestNeighbour)
        pszAlgName = "NearestNeighbour";
    else if (psWO->eResampleAlg == GRA_Bilinear)
        pszAlgName = "Bilinear";
    else if (psWO->eResampleAlg == GRA_Cubic)
        pszAlgName = "Cubic";
    else if (psWO->eResampleAlg == GRA_CubicSpline)
        pszAlgName = "CubicSpline";
    else if (psWO->eResampleAlg == GRA_Lanczos)
        pszAlgName = "Lanczos";
    else if (psWO->eResampleAlg == GRA_Average)
        pszAlgName = "Average";
    else if (psWO->eResampleAlg == GRA_RMS)
        pszAlgName = "RootMeanSquare";
    else if (psWO->eResampleAlg == GRA_Mode)
        pszAlgName = "Mode";
    else if (psWO->eResampleAlg == GRA_Max)
        pszAlgName = "Maximum";
    else if (psWO->eResampleAlg == GRA_Min)
        pszAlgName = "Minimum";
    else if (psWO->eResampleAlg == GRA_Med)
        pszAlgName = "Median";
    else if (psWO->eResampleAlg == GRA_Q1)
        pszAlgName = "Quartile1";
    else if (psWO->eResampleAlg == GRA_Q3)
        pszAlgName = "Quartile3";
    else if (psWO->eResampleAlg == GRA_Sum)
        pszAlgName = "Sum";
    else
        pszAlgName = "Unknown";

    CPLCreateXMLElementAndValue(psTree, "ResampleAlg", pszAlgName);

    /* -------------------------------------------------------------------- */
    /*      Working Data Type                                               */
    /* -------------------------------------------------------------------- */
    CPLCreateXMLElementAndValue(psTree, "WorkingDataType",
                                GDALGetDataTypeName(psWO->eWorkingDataType));

    /* -------------------------------------------------------------------- */
    /*      Name/value warp options.                                        */
    /* -------------------------------------------------------------------- */
    for (int iWO = 0; psWO->papszWarpOptions != nullptr &&
                      psWO->papszWarpOptions[iWO] != nullptr;
         iWO++)
    {
        char *pszName = nullptr;
        const char *pszValue =
            CPLParseNameValue(psWO->papszWarpOptions[iWO], &pszName);

        // EXTRA_ELTS is an internal detail that we will recover
        // no need to serialize it.
        // And CUTLINE is also serialized in a special way
        if (pszName != nullptr && !EQUAL(pszName, "EXTRA_ELTS") &&
            !EQUAL(pszName, "CUTLINE"))
        {
            CPLXMLNode *psOption =
                CPLCreateXMLElementAndValue(psTree, "Option", pszValue);

            CPLCreateXMLNode(CPLCreateXMLNode(psOption, CXT_Attribute, "name"),
                             CXT_Text, pszName);
        }

        CPLFree(pszName);
    }

    /* -------------------------------------------------------------------- */
    /*      Source and Destination Data Source                              */
    /* -------------------------------------------------------------------- */
    if (psWO->hSrcDS != nullptr)
    {
        CPLCreateXMLElementAndValue(psTree, "SourceDataset",
                                    GDALGetDescription(psWO->hSrcDS));

        CSLConstList papszOpenOptions =
            GDALDataset::FromHandle(psWO->hSrcDS)->GetOpenOptions();
        GDALSerializeOpenOptionsToXML(psTree, papszOpenOptions);
    }

    if (psWO->hDstDS != nullptr &&
        strlen(GDALGetDescription(psWO->hDstDS)) != 0)
    {
        CPLCreateXMLElementAndValue(psTree, "DestinationDataset",
                                    GDALGetDescription(psWO->hDstDS));
    }

    /* -------------------------------------------------------------------- */
    /*      Serialize transformer.                                          */
    /* -------------------------------------------------------------------- */
    if (psWO->pfnTransformer != nullptr)
    {
        CPLXMLNode *psTransformerContainer =
            CPLCreateXMLNode(psTree, CXT_Element, "Transformer");

        CPLXMLNode *psTransformerTree = GDALSerializeTransformer(
            psWO->pfnTransformer, psWO->pTransformerArg);

        if (psTransformerTree != nullptr)
            CPLAddXMLChild(psTransformerContainer, psTransformerTree);
    }

    /* -------------------------------------------------------------------- */
    /*      Band count and lists.                                           */
    /* -------------------------------------------------------------------- */
    CPLXMLNode *psBandList = nullptr;

    if (psWO->nBandCount != 0)
        psBandList = CPLCreateXMLNode(psTree, CXT_Element, "BandList");

    for (int i = 0; i < psWO->nBandCount; i++)
    {
        CPLXMLNode *psBand;

        psBand = CPLCreateXMLNode(psBandList, CXT_Element, "BandMapping");
        if (psWO->panSrcBands != nullptr)
            CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "src"),
                             CXT_Text,
                             CPLString().Printf("%d", psWO->panSrcBands[i]));
        if (psWO->panDstBands != nullptr)
            CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "dst"),
                             CXT_Text,
                             CPLString().Printf("%d", psWO->panDstBands[i]));

        if (psWO->padfSrcNoDataReal != nullptr)
        {
            CPLCreateXMLElementAndValue(
                psBand, "SrcNoDataReal",
                VRTSerializeNoData(psWO->padfSrcNoDataReal[i],
                                   psWO->eWorkingDataType, 16)
                    .c_str());
        }

        if (psWO->padfSrcNoDataImag != nullptr)
        {
            if (std::isnan(psWO->padfSrcNoDataImag[i]))
                CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "nan");
            else
                CPLCreateXMLElementAndValue(
                    psBand, "SrcNoDataImag",
                    CPLString().Printf("%.16g", psWO->padfSrcNoDataImag[i]));
        }
        // Compatibility with GDAL <= 2.2: if we serialize a SrcNoDataReal,
        // it needs a SrcNoDataImag as well
        else if (psWO->padfSrcNoDataReal != nullptr)
        {
            CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "0");
        }

        if (psWO->padfDstNoDataReal != nullptr)
        {
            CPLCreateXMLElementAndValue(
                psBand, "DstNoDataReal",
                VRTSerializeNoData(psWO->padfDstNoDataReal[i],
                                   psWO->eWorkingDataType, 16)
                    .c_str());
        }

        if (psWO->padfDstNoDataImag != nullptr)
        {
            if (std::isnan(psWO->padfDstNoDataImag[i]))
                CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "nan");
            else
                CPLCreateXMLElementAndValue(
                    psBand, "DstNoDataImag",
                    CPLString().Printf("%.16g", psWO->padfDstNoDataImag[i]));
        }
        // Compatibility with GDAL <= 2.2: if we serialize a DstNoDataReal,
        // it needs a SrcNoDataImag as well
        else if (psWO->padfDstNoDataReal != nullptr)
        {
            CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "0");
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Alpha bands.                                                    */
    /* -------------------------------------------------------------------- */
    if (psWO->nSrcAlphaBand > 0)
        CPLCreateXMLElementAndValue(
            psTree, "SrcAlphaBand",
            CPLString().Printf("%d", psWO->nSrcAlphaBand));

    if (psWO->nDstAlphaBand > 0)
        CPLCreateXMLElementAndValue(
            psTree, "DstAlphaBand",
            CPLString().Printf("%d", psWO->nDstAlphaBand));

    /* -------------------------------------------------------------------- */
    /*      Cutline.                                                        */
    /* -------------------------------------------------------------------- */
    if (psWO->hCutline != nullptr)
    {
        char *pszWKT = nullptr;
        if (OGR_G_ExportToWkt(static_cast<OGRGeometryH>(psWO->hCutline),
                              &pszWKT) == OGRERR_NONE)
        {
            CPLCreateXMLElementAndValue(psTree, "Cutline", pszWKT);
        }
        CPLFree(pszWKT);
    }

    if (psWO->dfCutlineBlendDist != 0.0)
        CPLCreateXMLElementAndValue(
            psTree, "CutlineBlendDist",
            CPLString().Printf("%.5g", psWO->dfCutlineBlendDist));

    return psTree;
}

/************************************************************************/
/*                     GDALDeserializeWarpOptions()                     */
/************************************************************************/

GDALWarpOptions *CPL_STDCALL GDALDeserializeWarpOptions(CPLXMLNode *psTree)

{
    CPLErrorReset();

    /* -------------------------------------------------------------------- */
    /*      Verify this is the right kind of object.                        */
    /* -------------------------------------------------------------------- */
    if (psTree == nullptr || psTree->eType != CXT_Element ||
        !EQUAL(psTree->pszValue, "GDALWarpOptions"))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Wrong node, unable to deserialize GDALWarpOptions.");
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create pre-initialized warp options.                            */
    /* -------------------------------------------------------------------- */
    GDALWarpOptions *psWO = GDALCreateWarpOptions();

    /* -------------------------------------------------------------------- */
    /*      Warp memory limit.                                              */
    /* -------------------------------------------------------------------- */
    psWO->dfWarpMemoryLimit =
        CPLAtof(CPLGetXMLValue(psTree, "WarpMemoryLimit", "0.0"));

    /* -------------------------------------------------------------------- */
    /*      resample algorithm                                              */
    /* -------------------------------------------------------------------- */
    const char *pszValue = CPLGetXMLValue(psTree, "ResampleAlg", "Default");

    if (EQUAL(pszValue, "NearestNeighbour"))
        psWO->eResampleAlg = GRA_NearestNeighbour;
    else if (EQUAL(pszValue, "Bilinear"))
        psWO->eResampleAlg = GRA_Bilinear;
    else if (EQUAL(pszValue, "Cubic"))
        psWO->eResampleAlg = GRA_Cubic;
    else if (EQUAL(pszValue, "CubicSpline"))
        psWO->eResampleAlg = GRA_CubicSpline;
    else if (EQUAL(pszValue, "Lanczos"))
        psWO->eResampleAlg = GRA_Lanczos;
    else if (EQUAL(pszValue, "Average"))
        psWO->eResampleAlg = GRA_Average;
    else if (EQUAL(pszValue, "RootMeanSquare"))
        psWO->eResampleAlg = GRA_RMS;
    else if (EQUAL(pszValue, "Mode"))
        psWO->eResampleAlg = GRA_Mode;
    else if (EQUAL(pszValue, "Maximum"))
        psWO->eResampleAlg = GRA_Max;
    else if (EQUAL(pszValue, "Minimum"))
        psWO->eResampleAlg = GRA_Min;
    else if (EQUAL(pszValue, "Median"))
        psWO->eResampleAlg = GRA_Med;
    else if (EQUAL(pszValue, "Quartile1"))
        psWO->eResampleAlg = GRA_Q1;
    else if (EQUAL(pszValue, "Quartile3"))
        psWO->eResampleAlg = GRA_Q3;
    else if (EQUAL(pszValue, "Sum"))
        psWO->eResampleAlg = GRA_Sum;
    else if (EQUAL(pszValue, "Default"))
        /* leave as is */;
    else
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Unrecognised ResampleAlg value '%s'.", pszValue);
    }

    /* -------------------------------------------------------------------- */
    /*      Working data type.                                              */
    /* -------------------------------------------------------------------- */
    psWO->eWorkingDataType = GDALGetDataTypeByName(
        CPLGetXMLValue(psTree, "WorkingDataType", "Unknown"));

    /* -------------------------------------------------------------------- */
    /*      Name/value warp options.                                        */
    /* -------------------------------------------------------------------- */
    for (CPLXMLNode *psItem = psTree->psChild; psItem != nullptr;
         psItem = psItem->psNext)
    {
        if (psItem->eType == CXT_Element && EQUAL(psItem->pszValue, "Option"))
        {
            const char *pszName = CPLGetXMLValue(psItem, "Name", nullptr);
            pszValue = CPLGetXMLValue(psItem, "", nullptr);

            if (pszName != nullptr && pszValue != nullptr)
            {
                psWO->papszWarpOptions =
                    CSLSetNameValue(psWO->papszWarpOptions, pszName, pszValue);
            }
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Source Dataset.                                                 */
    /* -------------------------------------------------------------------- */
    pszValue = CPLGetXMLValue(psTree, "SourceDataset", nullptr);

    if (pszValue != nullptr)
    {
        CPLXMLNode *psGeoLocNode =
            CPLSearchXMLNode(psTree, "GeoLocTransformer");
        if (psGeoLocNode)
        {
            CPLCreateXMLElementAndValue(psGeoLocNode, "SourceDataset",
                                        pszValue);
        }

        CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);

        char **papszOpenOptions = GDALDeserializeOpenOptionsFromXML(psTree);
        psWO->hSrcDS =
            GDALOpenEx(pszValue, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
                       nullptr, papszOpenOptions, nullptr);
        CSLDestroy(papszOpenOptions);
    }

    /* -------------------------------------------------------------------- */
    /*      Destination Dataset.                                            */
    /* -------------------------------------------------------------------- */
    pszValue = CPLGetXMLValue(psTree, "DestinationDataset", nullptr);

    if (pszValue != nullptr)
    {
        psWO->hDstDS = GDALOpenShared(pszValue, GA_Update);
    }

    /* -------------------------------------------------------------------- */
    /*      First, count band mappings so we can establish the bandcount.   */
    /* -------------------------------------------------------------------- */
    CPLXMLNode *psBandTree = CPLGetXMLNode(psTree, "BandList");

    int nBandCount = 0;
    CPLXMLNode *psBand = psBandTree ? psBandTree->psChild : nullptr;
    for (; psBand != nullptr; psBand = psBand->psNext)
    {
        if (psBand->eType != CXT_Element ||
            !EQUAL(psBand->pszValue, "BandMapping"))
            continue;

        nBandCount++;
    }

    GDALWarpInitDefaultBandMapping(psWO, nBandCount);

    /* ==================================================================== */
    /*      Now actually process each bandmapping.                          */
    /* ==================================================================== */
    int iBand = 0;

    psBand = psBandTree ? psBandTree->psChild : nullptr;

    for (; psBand != nullptr; psBand = psBand->psNext)
    {
        if (psBand->eType != CXT_Element ||
            !EQUAL(psBand->pszValue, "BandMapping"))
            continue;

        /* --------------------------------------------------------------------
         */
        /*      Source band */
        /* --------------------------------------------------------------------
         */
        pszValue = CPLGetXMLValue(psBand, "src", nullptr);
        if (pszValue != nullptr)
            psWO->panSrcBands[iBand] = atoi(pszValue);

        /* --------------------------------------------------------------------
         */
        /*      Destination band. */
        /* --------------------------------------------------------------------
         */
        pszValue = CPLGetXMLValue(psBand, "dst", nullptr);
        if (pszValue != nullptr)
            psWO->panDstBands[iBand] = atoi(pszValue);

        const auto NormalizeValue = [](const char *pszValueIn,
                                       GDALDataType eDataType) -> double
        {
            if (eDataType == GDT_Float32 &&
                CPLString().Printf("%.16g",
                                   static_cast<double>(
                                       std::numeric_limits<float>::lowest())) ==
                    pszValueIn)
            {
                return static_cast<double>(
                    std::numeric_limits<float>::lowest());
            }
            else if (eDataType == GDT_Float32 &&
                     CPLString().Printf(
                         "%.16g", static_cast<double>(
                                      std::numeric_limits<float>::max())) ==
                         pszValueIn)
            {
                return static_cast<double>(std::numeric_limits<float>::max());
            }
            else
            {
                return CPLAtof(pszValueIn);
            }
        };

        /* --------------------------------------------------------------------
         */
        /*      Source nodata. */
        /* --------------------------------------------------------------------
         */
        pszValue = CPLGetXMLValue(psBand, "SrcNoDataReal", nullptr);
        if (pszValue != nullptr)
        {
            GDALWarpInitSrcNoDataReal(psWO, -1.1e20);
            psWO->padfSrcNoDataReal[iBand] =
                NormalizeValue(pszValue, psWO->eWorkingDataType);
        }

        pszValue = CPLGetXMLValue(psBand, "SrcNoDataImag", nullptr);
        if (pszValue != nullptr)
        {
            GDALWarpInitSrcNoDataImag(psWO, 0);
            psWO->padfSrcNoDataImag[iBand] = CPLAtof(pszValue);
        }

        /* --------------------------------------------------------------------
         */
        /*      Destination nodata. */
        /* --------------------------------------------------------------------
         */
        pszValue = CPLGetXMLValue(psBand, "DstNoDataReal", nullptr);
        if (pszValue != nullptr)
        {
            GDALWarpInitDstNoDataReal(psWO, -1.1e20);
            psWO->padfDstNoDataReal[iBand] =
                NormalizeValue(pszValue, psWO->eWorkingDataType);
        }

        pszValue = CPLGetXMLValue(psBand, "DstNoDataImag", nullptr);
        if (pszValue != nullptr)
        {
            GDALWarpInitDstNoDataImag(psWO, 0);
            psWO->padfDstNoDataImag[iBand] = CPLAtof(pszValue);
        }

        iBand++;
    }

    /* -------------------------------------------------------------------- */
    /*      Alpha bands.                                                    */
    /* -------------------------------------------------------------------- */
    psWO->nSrcAlphaBand = atoi(CPLGetXMLValue(psTree, "SrcAlphaBand", "0"));
    psWO->nDstAlphaBand = atoi(CPLGetXMLValue(psTree, "DstAlphaBand", "0"));

    /* -------------------------------------------------------------------- */
    /*      Cutline.                                                        */
    /* -------------------------------------------------------------------- */
    const char *pszWKT = CPLGetXMLValue(psTree, "Cutline", nullptr);
    if (pszWKT)
    {
        char *pszWKTTemp = const_cast<char *>(pszWKT);
        OGRGeometryH hCutline = nullptr;
        OGR_G_CreateFromWkt(&pszWKTTemp, nullptr, &hCutline);
        psWO->hCutline = hCutline;
    }

    psWO->dfCutlineBlendDist =
        CPLAtof(CPLGetXMLValue(psTree, "CutlineBlendDist", "0"));

    /* -------------------------------------------------------------------- */
    /*      Transformation.                                                 */
    /* -------------------------------------------------------------------- */
    CPLXMLNode *psTransformer = CPLGetXMLNode(psTree, "Transformer");

    if (psTransformer != nullptr && psTransformer->psChild != nullptr)
    {
        GDALDeserializeTransformer(psTransformer->psChild,
                                   &(psWO->pfnTransformer),
                                   &(psWO->pTransformerArg));
    }

    /* -------------------------------------------------------------------- */
    /*      If any error has occurred, cleanup else return success.          */
    /* -------------------------------------------------------------------- */
    if (CPLGetLastErrorType() != CE_None)
    {
        if (psWO->pTransformerArg)
        {
            GDALDestroyTransformer(psWO->pTransformerArg);
            psWO->pTransformerArg = nullptr;
        }
        if (psWO->hSrcDS != nullptr)
        {
            GDALClose(psWO->hSrcDS);
            psWO->hSrcDS = nullptr;
        }
        if (psWO->hDstDS != nullptr)
        {
            GDALClose(psWO->hDstDS);
            psWO->hDstDS = nullptr;
        }
        GDALDestroyWarpOptions(psWO);
        return nullptr;
    }

    return psWO;
}

/************************************************************************/
/*                        GDALGetWarpResampleAlg()                      */
/************************************************************************/

/** Return a GDALResampleAlg from a string */
bool GDALGetWarpResampleAlg(const char *pszResampling,
                            GDALResampleAlg &eResampleAlg, bool bThrow)
{
    if (STARTS_WITH_CI(pszResampling, "near"))
        eResampleAlg = GRA_NearestNeighbour;
    else if (EQUAL(pszResampling, "bilinear"))
        eResampleAlg = GRA_Bilinear;
    else if (EQUAL(pszResampling, "cubic"))
        eResampleAlg = GRA_Cubic;
    else if (EQUAL(pszResampling, "cubicspline"))
        eResampleAlg = GRA_CubicSpline;
    else if (EQUAL(pszResampling, "lanczos"))
        eResampleAlg = GRA_Lanczos;
    else if (EQUAL(pszResampling, "average"))
        eResampleAlg = GRA_Average;
    else if (EQUAL(pszResampling, "rms"))
        eResampleAlg = GRA_RMS;
    else if (EQUAL(pszResampling, "mode"))
        eResampleAlg = GRA_Mode;
    else if (EQUAL(pszResampling, "max"))
        eResampleAlg = GRA_Max;
    else if (EQUAL(pszResampling, "min"))
        eResampleAlg = GRA_Min;
    else if (EQUAL(pszResampling, "med"))
        eResampleAlg = GRA_Med;
    else if (EQUAL(pszResampling, "q1"))
        eResampleAlg = GRA_Q1;
    else if (EQUAL(pszResampling, "q3"))
        eResampleAlg = GRA_Q3;
    else if (EQUAL(pszResampling, "sum"))
        eResampleAlg = GRA_Sum;
    else
    {
        if (bThrow)
        {
            throw std::invalid_argument("Unknown resampling method");
        }
        else
        {
            CPLError(CE_Failure, CPLE_IllegalArg,
                     "Unknown resampling method: %s.", pszResampling);
            return false;
        }
    }
    return true;
}