/******************************************************************************

 *

 * Project:  High Performance Image Reprojector

 * Purpose:  Implementation of the GDALWarpKernel class.  Implements the actual

 *           image warping for a "chunk" of input and output imagery already

 *           loaded into memory.

 * Author:   Frank Warmerdam, warmerdam@pobox.com

 *

 ******************************************************************************

 * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>

 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "cpl_port.h"

#include "gdalwarper.h"

#include <cfloat>

#include <cmath>

#include <cstddef>

#include <cstdlib>

#include <cstring>

#include <algorithm>

#include <limits>

#include <mutex>

#include <new>

#include <utility>

#include <vector>

#include "cpl_atomic_ops.h"

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_float.h"

#include "cpl_mask.h"

#include "cpl_multiproc.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "cpl_worker_thread_pool.h"

#include "cpl_quad_tree.h"

#include "gdal.h"

#include "gdal_alg.h"

#include "gdal_alg_priv.h"

#include "gdal_thread_pool.h"

#include "gdalresamplingkernels.h"

// #define CHECK_SUM_WITH_GEOS

#ifdef CHECK_SUM_WITH_GEOS

#include "ogr_geometry.h"

#include "ogr_geos.h"

#endif

#ifdef USE_NEON_OPTIMIZATIONS

#include "include_sse2neon.h"

#define USE_SSE2

#include "gdalsse_priv.h"

// We restrict to 64bit processors because they are guaranteed to have SSE2.

// Could possibly be used too on 32bit, but we would need to check at runtime.

#elif defined(__x86_64) || defined(_M_X64)

#define USE_SSE2

#include "gdalsse_priv.h"

#if __SSE4_1__

#include <smmintrin.h>

#endif

#if __SSE3__

#include <pmmintrin.h>

#endif

#endif

constexpr double BAND_DENSITY_THRESHOLD = 0.0000000001;

constexpr float SRC_DENSITY_THRESHOLD_FLOAT = 0.000000001f;

constexpr double SRC_DENSITY_THRESHOLD_DOUBLE = 0.000000001;

// #define INSTANTIATE_FLOAT64_SSE2_IMPL

static const int anGWKFilterRadius[] = {

    0,  // Nearest neighbour

    1,  // Bilinear

    2,  // Cubic Convolution (Catmull-Rom)

    2,  // Cubic B-Spline

    3,  // Lanczos windowed sinc

    0,  // Average

    0,  // Mode

    0,  // Reserved GRA_Gauss=7

    0,  // Max

    0,  // Min

    0,  // Med

    0,  // Q1

    0,  // Q3

    0,  // Sum

    0,  // RMS

};

static double GWKBilinear(double dfX);

static double GWKCubic(double dfX);

static double GWKBSpline(double dfX);

static double GWKLanczosSinc(double dfX);

static const FilterFuncType apfGWKFilter[] = {

    nullptr,         // Nearest neighbour

    GWKBilinear,     // Bilinear

    GWKCubic,        // Cubic Convolution (Catmull-Rom)

    GWKBSpline,      // Cubic B-Spline

    GWKLanczosSinc,  // Lanczos windowed sinc

    nullptr,         // Average

    nullptr,         // Mode

    nullptr,         // Reserved GRA_Gauss=7

    nullptr,         // Max

    nullptr,         // Min

    nullptr,         // Med

    nullptr,         // Q1

    nullptr,         // Q3

    nullptr,         // Sum

    nullptr,         // RMS

};

// TODO(schwehr): Can we make these functions have a const * const arg?

static double GWKBilinear4Values(double *padfVals);

static double GWKCubic4Values(double *padfVals);

static double GWKBSpline4Values(double *padfVals);

static double GWKLanczosSinc4Values(double *padfVals);

static const FilterFunc4ValuesType apfGWKFilter4Values[] = {

    nullptr,                // Nearest neighbour

    GWKBilinear4Values,     // Bilinear

    GWKCubic4Values,        // Cubic Convolution (Catmull-Rom)

    GWKBSpline4Values,      // Cubic B-Spline

    GWKLanczosSinc4Values,  // Lanczos windowed sinc

    nullptr,                // Average

    nullptr,                // Mode

    nullptr,                // Reserved GRA_Gauss=7

    nullptr,                // Max

    nullptr,                // Min

    nullptr,                // Med

    nullptr,                // Q1

    nullptr,                // Q3

    nullptr,                // Sum

    nullptr,                // RMS

};

int GWKGetFilterRadius(GDALResampleAlg eResampleAlg)

{

    static_assert(CPL_ARRAYSIZE(anGWKFilterRadius) == GRA_LAST_VALUE + 1,

                  "Bad size of anGWKFilterRadius");

    return anGWKFilterRadius[eResampleAlg];

}

FilterFuncType GWKGetFilterFunc(GDALResampleAlg eResampleAlg)

{

    static_assert(CPL_ARRAYSIZE(apfGWKFilter) == GRA_LAST_VALUE + 1,

                  "Bad size of apfGWKFilter");

    return apfGWKFilter[eResampleAlg];

}

FilterFunc4ValuesType GWKGetFilterFunc4Values(GDALResampleAlg eResampleAlg)

{

    static_assert(CPL_ARRAYSIZE(apfGWKFilter4Values) == GRA_LAST_VALUE + 1,

                  "Bad size of apfGWKFilter4Values");

    return apfGWKFilter4Values[eResampleAlg];

}

static CPLErr GWKGeneralCase(GDALWarpKernel *);

static CPLErr GWKRealCase(GDALWarpKernel *poWK);

static CPLErr GWKNearestNoMasksOrDstDensityOnlyByte(GDALWarpKernel *poWK);

static CPLErr GWKBilinearNoMasksOrDstDensityOnlyByte(GDALWarpKernel *poWK);

static CPLErr GWKCubicNoMasksOrDstDensityOnlyByte(GDALWarpKernel *poWK);

static CPLErr GWKCubicNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK);

#ifdef INSTANTIATE_FLOAT64_SSE2_IMPL

static CPLErr GWKCubicNoMasksOrDstDensityOnlyDouble(GDALWarpKernel *poWK);

#endif

static CPLErr GWKCubicSplineNoMasksOrDstDensityOnlyByte(GDALWarpKernel *poWK);

static CPLErr GWKNearestByte(GDALWarpKernel *poWK);

static CPLErr GWKNearestNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);

static CPLErr GWKBilinearNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);

static CPLErr GWKBilinearNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK);

#ifdef INSTANTIATE_FLOAT64_SSE2_IMPL

static CPLErr GWKBilinearNoMasksOrDstDensityOnlyDouble(GDALWarpKernel *poWK);

#endif

static CPLErr GWKCubicNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);

static CPLErr GWKCubicSplineNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);

static CPLErr GWKNearestShort(GDALWarpKernel *poWK);

static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK);

static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK);

static CPLErr GWKNearestFloat(GDALWarpKernel *poWK);

static CPLErr GWKAverageOrMode(GDALWarpKernel *);

static CPLErr GWKSumPreserving(GDALWarpKernel *);

static CPLErr GWKCubicNoMasksOrDstDensityOnlyUShort(GDALWarpKernel *);

static CPLErr GWKCubicSplineNoMasksOrDstDensityOnlyUShort(GDALWarpKernel *);

static CPLErr GWKBilinearNoMasksOrDstDensityOnlyUShort(GDALWarpKernel *);

/************************************************************************/

/*                             GWKJobStruct                             */

/************************************************************************/

struct GWKJobStruct

{

    std::mutex &mutex;

    std::condition_variable &cv;

    int counterSingleThreaded = 0;

    int &counter;

    bool &stopFlag;

    GDALWarpKernel *poWK = nullptr;

    int iYMin = 0;

    int iYMax = 0;

    int (*pfnProgress)(GWKJobStruct *psJob) = nullptr;

    void *pTransformerArg = nullptr;

    // used by GWKRun() to assign the proper pTransformerArg

    void (*pfnFunc)(void *) = nullptr;

    GWKJobStruct(std::mutex &mutex_, std::condition_variable &cv_,

                 int &counter_, bool &stopFlag_)

        : mutex(mutex_), cv(cv_), counter(counter_), stopFlag(stopFlag_)

    {

    }

};

struct GWKThreadData

{

    std::unique_ptr<CPLJobQueue> poJobQueue{};

    std::unique_ptr<std::vector<GWKJobStruct>> threadJobs{};

    int nMaxThreads{0};

    int counter{0};

    bool stopFlag{false};

    std::mutex mutex{};

    std::condition_variable cv{};

    bool bTransformerArgInputAssignedToThread{false};

    void *pTransformerArgInput{

        nullptr};  // owned by calling layer. Not to be destroyed

    std::map<GIntBig, void *> mapThreadToTransformerArg{};

    int nTotalThreadCountForThisRun = 0;

    int nCurThreadCountForThisRun = 0;

};

/************************************************************************/

/*                         GWKProgressThread()                          */

/************************************************************************/

// Return TRUE if the computation must be interrupted.

static int GWKProgressThread(GWKJobStruct *psJob)

{

    bool stop = false;

    {

        std::lock_guard<std::mutex> lock(psJob->mutex);

        psJob->counter++;

        stop = psJob->stopFlag;

    }

    psJob->cv.notify_one();

    return stop;

}

/************************************************************************/

/*                       GWKProgressMonoThread()                        */

/************************************************************************/

// Return TRUE if the computation must be interrupted.

static int GWKProgressMonoThread(GWKJobStruct *psJob)

{

    GDALWarpKernel *poWK = psJob->poWK;

    if (!poWK->pfnProgress(poWK->dfProgressBase +

                               poWK->dfProgressScale *

                                   (++psJob->counterSingleThreaded /

                                    static_cast<double>(psJob->iYMax)),

                           "", poWK->pProgress))

    {

        CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

        psJob->stopFlag = true;

        return TRUE;

    }

    return FALSE;

}

/************************************************************************/

/*                        GWKGenericMonoThread()                        */

/************************************************************************/

static CPLErr GWKGenericMonoThread(GDALWarpKernel *poWK,

                                   void (*pfnFunc)(void *pUserData))

{

    GWKThreadData td;

    // NOTE: the mutex is not used.

    GWKJobStruct job(td.mutex, td.cv, td.counter, td.stopFlag);

    job.poWK = poWK;

    job.iYMin = 0;

    job.iYMax = poWK->nDstYSize;

    job.pfnProgress = GWKProgressMonoThread;

    job.pTransformerArg = poWK->pTransformerArg;

    job.counterSingleThreaded = td.counter;

    pfnFunc(&job);

    td.counter = job.counterSingleThreaded;

    return td.stopFlag ? CE_Failure : CE_None;

}

/************************************************************************/

/*                          GWKThreadsCreate()                          */

/************************************************************************/

void *GWKThreadsCreate(char **papszWarpOptions,

                       GDALTransformerFunc /* pfnTransformer */,

                       void *pTransformerArg)

{

    const int nThreads = GDALGetNumThreads(papszWarpOptions, "NUM_THREADS",

                                           GDAL_DEFAULT_MAX_THREAD_COUNT,

                                           /* bDefaultAllCPUs = */ false);

    GWKThreadData *psThreadData = new GWKThreadData();

    auto poThreadPool =

        nThreads > 1 ? GDALGetGlobalThreadPool(nThreads) : nullptr;

    if (poThreadPool)

    {

        psThreadData->nMaxThreads = nThreads;

        psThreadData->threadJobs.reset(new std::vector<GWKJobStruct>(

            nThreads,

            GWKJobStruct(psThreadData->mutex, psThreadData->cv,

                         psThreadData->counter, psThreadData->stopFlag)));

        psThreadData->poJobQueue = poThreadPool->CreateJobQueue();

        psThreadData->pTransformerArgInput = pTransformerArg;

    }

    return psThreadData;

}

/************************************************************************/

/*                           GWKThreadsEnd()                            */

/************************************************************************/

void GWKThreadsEnd(void *psThreadDataIn)

{

    if (psThreadDataIn == nullptr)

        return;

    GWKThreadData *psThreadData = static_cast<GWKThreadData *>(psThreadDataIn);

    if (psThreadData->poJobQueue)

    {

        // cppcheck-suppress constVariableReference

        for (auto &pair : psThreadData->mapThreadToTransformerArg)

        {

            CPLAssert(pair.second != psThreadData->pTransformerArgInput);

            GDALDestroyTransformer(pair.second);

        }

        psThreadData->poJobQueue.reset();

    }

    delete psThreadData;

}

/************************************************************************/

/*                         ThreadFuncAdapter()                          */

/************************************************************************/

static void ThreadFuncAdapter(void *pData)

{

    GWKJobStruct *psJob = static_cast<GWKJobStruct *>(pData);

    GWKThreadData *psThreadData =

        static_cast<GWKThreadData *>(psJob->poWK->psThreadData);

    // Look if we have already a per-thread transformer

    void *pTransformerArg = nullptr;

    const GIntBig nThreadId = CPLGetPID();

    {

        std::lock_guard<std::mutex> lock(psThreadData->mutex);

        ++psThreadData->nCurThreadCountForThisRun;

        auto oIter = psThreadData->mapThreadToTransformerArg.find(nThreadId);

        if (oIter != psThreadData->mapThreadToTransformerArg.end())

        {

            pTransformerArg = oIter->second;

        }

        else if (!psThreadData->bTransformerArgInputAssignedToThread &&

                 psThreadData->nCurThreadCountForThisRun ==

                     psThreadData->nTotalThreadCountForThisRun)

        {

            // If we are the last thread to be started, temporarily borrow the

            // original transformer

            psThreadData->bTransformerArgInputAssignedToThread = true;

            pTransformerArg = psThreadData->pTransformerArgInput;

            psThreadData->mapThreadToTransformerArg[nThreadId] =

                pTransformerArg;

        }

        if (pTransformerArg == nullptr)

        {

            CPLAssert(psThreadData->pTransformerArgInput != nullptr);

            CPLAssert(!psThreadData->bTransformerArgInputAssignedToThread);

        }

    }

    // If no transformer assigned to current thread, instantiate one

    if (pTransformerArg == nullptr)

    {

        // This somehow assumes that GDALCloneTransformer() is thread-safe

        // which should normally be the case.

        pTransformerArg =

            GDALCloneTransformer(psThreadData->pTransformerArgInput);

        // Lock for the stop flag and the transformer map.

        std::lock_guard<std::mutex> lock(psThreadData->mutex);

        if (!pTransformerArg)

        {

            psJob->stopFlag = true;

            return;

        }

        psThreadData->mapThreadToTransformerArg[nThreadId] = pTransformerArg;

    }

    psJob->pTransformerArg = pTransformerArg;

    psJob->pfnFunc(pData);

    // Give back original transformer, if borrowed.

    {

        std::lock_guard<std::mutex> lock(psThreadData->mutex);

        if (psThreadData->bTransformerArgInputAssignedToThread &&

            pTransformerArg == psThreadData->pTransformerArgInput)

        {

            psThreadData->mapThreadToTransformerArg.erase(

                psThreadData->mapThreadToTransformerArg.find(nThreadId));

            psThreadData->bTransformerArgInputAssignedToThread = false;

        }

    }

}

/************************************************************************/

/*                               GWKRun()                               */

/************************************************************************/

static CPLErr GWKRun(GDALWarpKernel *poWK, const char *pszFuncName,

                     void (*pfnFunc)(void *pUserData))

{

    const int nDstYSize = poWK->nDstYSize;

    CPLDebug("GDAL",

             "GDALWarpKernel()::%s() "

             "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",

             pszFuncName, poWK->nSrcXOff, poWK->nSrcYOff, poWK->nSrcXSize,

             poWK->nSrcYSize, poWK->nDstXOff, poWK->nDstYOff, poWK->nDstXSize,

             poWK->nDstYSize);

    if (!poWK->pfnProgress(poWK->dfProgressBase, "", poWK->pProgress))

    {

        CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

        return CE_Failure;

    }

    GWKThreadData *psThreadData =

        static_cast<GWKThreadData *>(poWK->psThreadData);

    if (psThreadData == nullptr || psThreadData->poJobQueue == nullptr)

    {

        return GWKGenericMonoThread(poWK, pfnFunc);

    }

    int nThreads = std::min(psThreadData->nMaxThreads, nDstYSize / 2);

    // Config option mostly useful for tests to be able to test multithreading

    // with small rasters

    const int nWarpChunkSize =

        atoi(CPLGetConfigOption("WARP_THREAD_CHUNK_SIZE", "65536"));

    if (nWarpChunkSize > 0)

    {

        GIntBig nChunks =

            static_cast<GIntBig>(nDstYSize) * poWK->nDstXSize / nWarpChunkSize;

        if (nThreads > nChunks)

            nThreads = static_cast<int>(nChunks);

    }

    if (nThreads <= 0)

        nThreads = 1;

    CPLDebug("WARP", "Using %d threads", nThreads);

    auto &jobs = *psThreadData->threadJobs;

    CPLAssert(static_cast<int>(jobs.size()) >= nThreads);

    // Fill-in job structures.

    for (int i = 0; i < nThreads; ++i)

    {

        auto &job = jobs[i];

        job.poWK = poWK;

        job.iYMin =

            static_cast<int>(static_cast<int64_t>(i) * nDstYSize / nThreads);

        job.iYMax = static_cast<int>(static_cast<int64_t>(i + 1) * nDstYSize /

                                     nThreads);

        if (poWK->pfnProgress != GDALDummyProgress)

            job.pfnProgress = GWKProgressThread;

        job.pfnFunc = pfnFunc;

    }

    bool bStopFlag;

    {

        std::unique_lock<std::mutex> lock(psThreadData->mutex);

        psThreadData->nTotalThreadCountForThisRun = nThreads;

        // coverity[missing_lock]

        psThreadData->nCurThreadCountForThisRun = 0;

        // Start jobs.

        for (int i = 0; i < nThreads; ++i)

        {

            auto &job = jobs[i];

            psThreadData->poJobQueue->SubmitJob(ThreadFuncAdapter,

                                                static_cast<void *>(&job));

        }

        /* --------------------------------------------------------------------

         */

        /*      Report progress. */

        /* --------------------------------------------------------------------

         */

        if (poWK->pfnProgress != GDALDummyProgress)

        {

            while (psThreadData->counter < nDstYSize)

            {

                psThreadData->cv.wait(lock);

                if (!poWK->pfnProgress(poWK->dfProgressBase +

                                           poWK->dfProgressScale *

                                               (psThreadData->counter /

                                                static_cast<double>(nDstYSize)),

                                       "", poWK->pProgress))

                {

                    CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

                    psThreadData->stopFlag = true;

                    break;

                }

            }

        }

        bStopFlag = psThreadData->stopFlag;

    }

    /* -------------------------------------------------------------------- */

    /*      Wait for all jobs to complete.                                  */

    /* -------------------------------------------------------------------- */

    psThreadData->poJobQueue->WaitCompletion();

    return bStopFlag ? CE_Failure : CE_None;

}

/************************************************************************/

/* ==================================================================== */

/*                            GDALWarpKernel                            */

/* ==================================================================== */

/************************************************************************/

/**

 * \class GDALWarpKernel "gdalwarper.h"

 *

 * Low level image warping class.

 *

 * This class is responsible for low level image warping for one

 * "chunk" of imagery.  The class is essentially a structure with all

 * data members public - primarily so that new special-case functions

 * can be added without changing the class declaration.

 *

 * Applications are normally intended to interactive with warping facilities

 * through the GDALWarpOperation class, though the GDALWarpKernel can in

 * theory be used directly if great care is taken in setting up the

 * control data.

 *

 * <h3>Design Issues</h3>

 *

 * The intention is that PerformWarp() would analyze the setup in terms

 * of the datatype, resampling type, and validity/density mask usage and

 * pick one of many specific implementations of the warping algorithm over

 * a continuum of optimization vs. generality.  At one end there will be a

 * reference general purpose implementation of the algorithm that supports

 * any data type (working internally in double precision complex), all three

 * resampling types, and any or all of the validity/density masks.  At the

 * other end would be highly optimized algorithms for common cases like

 * nearest neighbour resampling on GDT_UInt8 data with no masks.

 *

 * The full set of optimized versions have not been decided but we should

 * expect to have at least:

 *  - One for each resampling algorithm for 8bit data with no masks.

 *  - One for each resampling algorithm for float data with no masks.

 *  - One for each resampling algorithm for float data with any/all masks

 *    (essentially the generic case for just float data).

 *  - One for each resampling algorithm for 8bit data with support for

 *    input validity masks (per band or per pixel).  This handles the common

 *    case of nodata masking.

 *  - One for each resampling algorithm for float data with support for

 *    input validity masks (per band or per pixel).  This handles the common

 *    case of nodata masking.

 *

 * Some of the specializations would operate on all bands in one pass

 * (especially the ones without masking would do this), while others might

 * process each band individually to reduce code complexity.

 *

 * <h3>Masking Semantics</h3>

 *

 * A detailed explanation of the semantics of the validity and density masks,

 * and their effects on resampling kernels is needed here.

 */

/************************************************************************/

/*                     GDALWarpKernel Data Members                      */

/************************************************************************/

/**

 * \var GDALResampleAlg GDALWarpKernel::eResample;

 *

 * Resampling algorithm.

 *

 * The resampling algorithm to use.  One of GRA_NearestNeighbour, GRA_Bilinear,

 * GRA_Cubic, GRA_CubicSpline, GRA_Lanczos, GRA_Average, GRA_RMS,

 * GRA_Mode or GRA_Sum.

 *

 * This field is required. GDT_NearestNeighbour may be used as a default

 * value.

 */

/**

 * \var GDALDataType GDALWarpKernel::eWorkingDataType;

 *

 * Working pixel data type.

 *

 * The datatype of pixels in the source image (papabySrcimage) and

 * destination image (papabyDstImage) buffers.  Note that operations on

 * some data types (such as GDT_UInt8) may be much better optimized than other

 * less common cases.

 *

 * This field is required.  It may not be GDT_Unknown.

 */

/**

 * \var int GDALWarpKernel::nBands;

 *

 * Number of bands.

 *

 * The number of bands (layers) of imagery being warped.  Determines the

 * number of entries in the papabySrcImage, papanBandSrcValid,

 * and papabyDstImage arrays.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nSrcXSize;

 *

 * Source image width in pixels.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nSrcYSize;

 *

 * Source image height in pixels.

 *

 * This field is required.

 */

/**

 * \var double GDALWarpKernel::dfSrcXExtraSize;

 *

 * Number of pixels included in nSrcXSize that are present on the edges of

 * the area of interest to take into account the width of the kernel.

 *

 * This field is required.

 */

/**

 * \var double GDALWarpKernel::dfSrcYExtraSize;

 *

 * Number of pixels included in nSrcYExtraSize that are present on the edges of

 * the area of interest to take into account the height of the kernel.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::papabySrcImage;

 *

 * Array of source image band data.

 *

 * This is an array of pointers (of size GDALWarpKernel::nBands) pointers

 * to image data.  Each individual band of image data is organized as a single

 * block of image data in left to right, then bottom to top order.  The actual

 * type of the image data is determined by GDALWarpKernel::eWorkingDataType.

 *

 * To access the pixel value for the (x=3, y=4) pixel (zero based) of

 * the second band with eWorkingDataType set to GDT_Float32 use code like

 * this:

 *

 * \code

 *   float dfPixelValue;

 *   int   nBand = 2-1;  // Band indexes are zero based.

 *   int   nPixel = 3; // Zero based.

 *   int   nLine = 4;  // Zero based.

 *

 *   assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );

 *   assert( nLine >= 0 && nLine < poKern->nSrcYSize );

 *   assert( nBand >= 0 && nBand < poKern->nBands );

 *   dfPixelValue = ((float *) poKern->papabySrcImage[nBand])

 *                                  [nPixel + nLine * poKern->nSrcXSize];

 * \endcode

 *

 * This field is required.

 */

/**

 * \var GUInt32 **GDALWarpKernel::papanBandSrcValid;

 *

 * Per band validity mask for source pixels.

 *

 * Array of pixel validity mask layers for each source band.   Each of

 * the mask layers is the same size (in pixels) as the source image with

 * one bit per pixel.  Note that it is legal (and common) for this to be

 * NULL indicating that none of the pixels are invalidated, or for some

 * band validity masks to be NULL in which case all pixels of the band are

 * valid.  The following code can be used to test the validity of a particular

 * pixel.

 *

 * \code

 *   int   bIsValid = TRUE;

 *   int   nBand = 2-1;  // Band indexes are zero based.

 *   int   nPixel = 3; // Zero based.

 *   int   nLine = 4;  // Zero based.

 *

 *   assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );

 *   assert( nLine >= 0 && nLine < poKern->nSrcYSize );

 *   assert( nBand >= 0 && nBand < poKern->nBands );

 *

 *   if( poKern->papanBandSrcValid != NULL

 *       && poKern->papanBandSrcValid[nBand] != NULL )

 *   {

 *       GUInt32 *panBandMask = poKern->papanBandSrcValid[nBand];

 *       int    iPixelOffset = nPixel + nLine * poKern->nSrcXSize;

 *

 *       bIsValid = CPLMaskGet(panBandMask, iPixelOffset)

 *   }

 * \endcode

 */

/**

 * \var GUInt32 *GDALWarpKernel::panUnifiedSrcValid;

 *

 * Per pixel validity mask for source pixels.

 *

 * A single validity mask layer that applies to the pixels of all source

 * bands.  It is accessed similarly to papanBandSrcValid, but without the

 * extra level of band indirection.

 *

 * This pointer may be NULL indicating that all pixels are valid.

 *

 * Note that if both panUnifiedSrcValid, and papanBandSrcValid are available,

 * the pixel isn't considered to be valid unless both arrays indicate it is

 * valid.

 */

/**

 * \var float *GDALWarpKernel::pafUnifiedSrcDensity;

 *

 * Per pixel density mask for source pixels.

 *

 * A single density mask layer that applies to the pixels of all source

 * bands.  It contains values between 0.0 and 1.0 indicating the degree to

 * which this pixel should be allowed to contribute to the output result.

 *

 * This pointer may be NULL indicating that all pixels have a density of 1.0.

 *

 * The density for a pixel may be accessed like this:

 *

 * \code

 *   float fDensity = 1.0;

 *   int nPixel = 3;  // Zero based.

 *   int nLine = 4;   // Zero based.

 *

 *   assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );

 *   assert( nLine >= 0 && nLine < poKern->nSrcYSize );

 *   if( poKern->pafUnifiedSrcDensity != NULL )

 *     fDensity = poKern->pafUnifiedSrcDensity

 *                                  [nPixel + nLine * poKern->nSrcXSize];

 * \endcode

 */

/**

 * \var int GDALWarpKernel::nDstXSize;

 *

 * Width of destination image in pixels.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nDstYSize;

 *

 * Height of destination image in pixels.

 *

 * This field is required.

 */

/**

 * \var GByte **GDALWarpKernel::papabyDstImage;

 *

 * Array of destination image band data.

 *

 * This is an array of pointers (of size GDALWarpKernel::nBands) pointers

 * to image data.  Each individual band of image data is organized as a single

 * block of image data in left to right, then bottom to top order.  The actual

 * type of the image data is determined by GDALWarpKernel::eWorkingDataType.

 *

 * To access the pixel value for the (x=3, y=4) pixel (zero based) of

 * the second band with eWorkingDataType set to GDT_Float32 use code like

 * this:

 *

 * \code

 *   float dfPixelValue;

 *   int   nBand = 2-1;  // Band indexes are zero based.

 *   int   nPixel = 3; // Zero based.

 *   int   nLine = 4;  // Zero based.

 *

 *   assert( nPixel >= 0 && nPixel < poKern->nDstXSize );

 *   assert( nLine >= 0 && nLine < poKern->nDstYSize );

 *   assert( nBand >= 0 && nBand < poKern->nBands );

 *   dfPixelValue = ((float *) poKern->papabyDstImage[nBand])

 *                                  [nPixel + nLine * poKern->nSrcYSize];

 * \endcode

 *

 * This field is required.

 */

/**

 * \var GUInt32 *GDALWarpKernel::panDstValid;

 *

 * Per pixel validity mask for destination pixels.

 *

 * A single validity mask layer that applies to the pixels of all destination

 * bands.  It is accessed similarly to papanUnitifiedSrcValid, but based

 * on the size of the destination image.

 *

 * This pointer may be NULL indicating that all pixels are valid.

 */

/**

 * \var float *GDALWarpKernel::pafDstDensity;

 *

 * Per pixel density mask for destination pixels.

 *

 * A single density mask layer that applies to the pixels of all destination

 * bands.  It contains values between 0.0 and 1.0.

 *

 * This pointer may be NULL indicating that all pixels have a density of 1.0.

 *

 * The density for a pixel may be accessed like this:

 *

 * \code

 *   float fDensity = 1.0;

 *   int   nPixel = 3; // Zero based.

 *   int   nLine = 4;  // Zero based.

 *

 *   assert( nPixel >= 0 && nPixel < poKern->nDstXSize );

 *   assert( nLine >= 0 && nLine < poKern->nDstYSize );

 *   if( poKern->pafDstDensity != NULL )

 *     fDensity = poKern->pafDstDensity[nPixel + nLine * poKern->nDstXSize];

 * \endcode

 */

/**

 * \var int GDALWarpKernel::nSrcXOff;

 *

 * X offset to source pixel coordinates for transformation.

 *

 * See pfnTransformer.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nSrcYOff;

 *

 * Y offset to source pixel coordinates for transformation.

 *

 * See pfnTransformer.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nDstXOff;

 *

 * X offset to destination pixel coordinates for transformation.

 *

 * See pfnTransformer.

 *

 * This field is required.

 */

/**

 * \var int GDALWarpKernel::nDstYOff;

 *

 * Y offset to destination pixel coordinates for transformation.

 *

 * See pfnTransformer.

 *

 * This field is required.

 */

/**

 * \var GDALTransformerFunc GDALWarpKernel::pfnTransformer;

 *

 * Source/destination location transformer.

 *

 * The function to call to transform coordinates between source image

 * pixel/line coordinates and destination image pixel/line coordinates.

 * See GDALTransformerFunc() for details of the semantics of this function.

 *

 * The GDALWarpKern algorithm will only ever use this transformer in

 * "destination to source" mode (bDstToSrc=TRUE), and will always pass

 * partial or complete scanlines of points in the destination image as

 * input.  This means, among other things, that it is safe to the

 * approximating transform GDALApproxTransform() as the transformation

 * function.

 *

 * Source and destination images may be subsets of a larger overall image.

 * The transformation algorithms will expect and return pixel/line coordinates

 * in terms of this larger image, so coordinates need to be offset by

 * the offsets specified in nSrcXOff, nSrcYOff, nDstXOff, and nDstYOff before

 * passing to pfnTransformer, and after return from it.

 *

 * The GDALWarpKernel::pfnTransformerArg value will be passed as the callback

 * data to this function when it is called.

 *

 * This field is required.

 */

/**

 * \var void *GDALWarpKernel::pTransformerArg;

 *

 * Callback data for pfnTransformer.

 *

 * This field may be NULL if not required for the pfnTransformer being used.

 */

/**

 * \var GDALProgressFunc GDALWarpKernel::pfnProgress;

 *

 * The function to call to report progress of the algorithm, and to check

 * for a requested termination of the operation.  It operates according to

 * GDALProgressFunc() semantics.

 *

 * Generally speaking the progress function will be invoked for each

 * scanline of the destination buffer that has been processed.

 *

 * This field may be NULL (internally set to GDALDummyProgress()).

 */

/**

 * \var void *GDALWarpKernel::pProgress;

 *

 * Callback data for pfnProgress.

 *

 * This field may be NULL if not required for the pfnProgress being used.

 */

/************************************************************************/

/*                           GDALWarpKernel()                           */

/************************************************************************/

GDALWarpKernel::GDALWarpKernel()

    : papszWarpOptions(nullptr), eResample(GRA_NearestNeighbour),

      eWorkingDataType(GDT_Unknown), nBands(0), nSrcXSize(0), nSrcYSize(0),

      dfSrcXExtraSize(0.0), dfSrcYExtraSize(0.0), papabySrcImage(nullptr),

      papanBandSrcValid(nullptr), panUnifiedSrcValid(nullptr),

      pafUnifiedSrcDensity(nullptr), nDstXSize(0), nDstYSize(0),

      papabyDstImage(nullptr), panDstValid(nullptr), pafDstDensity(nullptr),

      dfXScale(1.0), dfYScale(1.0), dfXFilter(0.0), dfYFilter(0.0), nXRadius(0),

      nYRadius(0), nFiltInitX(0), nFiltInitY(0), nSrcXOff(0), nSrcYOff(0),

      nDstXOff(0), nDstYOff(0), pfnTransformer(nullptr),

      pTransformerArg(nullptr), pfnProgress(GDALDummyProgress),

      pProgress(nullptr), dfProgressBase(0.0), dfProgressScale(1.0),

      padfDstNoDataReal(nullptr), psThreadData(nullptr),

      eTieStrategy(GWKTS_First)

{

}

/************************************************************************/

/*                          ~GDALWarpKernel()                           */

/************************************************************************/

GDALWarpKernel::~GDALWarpKernel()

{

}

/************************************************************************/

/*                              getArea()                               */

/************************************************************************/

typedef std::pair<double, double> XYPair;

typedef std::vector<XYPair> XYPoly;

// poly may or may not be closed.

static double getArea(const XYPoly &poly)

{

    // CPLAssert(poly.size() >= 2);

    const size_t nPointCount = poly.size();

    double dfAreaSum =

        poly[0].first * (poly[1].second - poly[nPointCount - 1].second);

    for (size_t i = 1; i < nPointCount - 1; i++)

    {

        dfAreaSum += poly[i].first * (poly[i + 1].second - poly[i - 1].second);

    }

    dfAreaSum += poly[nPointCount - 1].first *

                 (poly[0].second - poly[nPointCount - 2].second);

    return 0.5 * std::fabs(dfAreaSum);

}

/************************************************************************/

/*                       CanUse4SamplesFormula()                        */

/************************************************************************/

static bool CanUse4SamplesFormula(const GDALWarpKernel *poWK)

{

    if (poWK->eResample == GRA_Bilinear || poWK->eResample == GRA_Cubic)

    {

        // Use 4-sample formula if we are not downsampling by more than a

        // factor of 1:2

        if (poWK->dfXScale > 0.5 && poWK->dfYScale > 0.5)

            return true;

        CPLDebugOnce("WARP",

                     "Not using 4-sample bilinear/bicubic formula because "

                     "XSCALE(=%f) and/or YSCALE(=%f) <= 0.5",

                     poWK->dfXScale, poWK->dfYScale);

    }

    return false;

}