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

 *

 * Project:  GDAL Utilities

 * Purpose:  Footprint OGR shapes into a GDAL raster.

 * Authors:  Even Rouault, <even dot rouault at spatialys dot com>

 *

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

 * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot com>

 *

 * SPDX-License-Identifier: MIT

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

#include "cpl_port.h"

#include "gdal_utils.h"

#include "gdal_utils_priv.h"

#include <cmath>

#include <cstdio>

#include <cstdlib>

#include <cstring>

#include <algorithm>

#include <limits>

#include <vector>

#include "commonutils.h"

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "gdal.h"

#include "gdal_alg.h"

#include "gdal_priv.h"

#include "gdal_maskbands.h"

#include "ogr_api.h"

#include "ogr_core.h"

#include "memdataset.h"

#include "ogrsf_frmts.h"

#include "ogr_spatialref.h"

#include "gdalargumentparser.h"

constexpr const char *DEFAULT_LAYER_NAME = "footprint";

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

/*                         GDALFootprintOptions                         */

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

struct GDALFootprintOptions

{

    /*! output format. Use the short format name. */

    std::string osFormat{};

    /*! the progress function to use */

    GDALProgressFunc pfnProgress = GDALDummyProgress;

    /*! pointer to the progress data variable */

    void *pProgressData = nullptr;

    bool bCreateOutput = false;

    std::string osDestLayerName{};

    /*! Layer creation options */

    CPLStringList aosLCO{};

    /*! Dataset creation options */

    CPLStringList aosDSCO{};

    /*! Overview index: 0 = first overview level */

    int nOvrIndex = -1;

    /** Whether output geometry should be in georeferenced coordinates, if

     * possible (if explicitly requested, bOutCSGeorefRequested is also set)

     * false = in pixel coordinates

     */

    bool bOutCSGeoref = true;

    /** Whether -t_cs georef has been explicitly set */

    bool bOutCSGeorefRequested = false;

    OGRSpatialReference oOutputSRS{};

    bool bSplitPolys = false;

    double dfDensifyDistance = 0;

    double dfSimplifyTolerance = 0;

    bool bConvexHull = false;

    double dfMinRingArea = 0;

    int nMaxPoints = 100;

    /*! Source bands to take into account */

    std::vector<int> anBands{};

    /*! Whether to combine bands unioning (true) or intersecting (false) */

    bool bCombineBandsUnion = true;

    /*! Field name where to write the path of the raster. Empty if not desired */

    std::string osLocationFieldName = "location";

    /*! Clears the osLocationFieldName var when set */

    bool bClearLocation = false;

    /*! Whether to force writing absolute paths in location field. */

    bool bAbsolutePath = false;

    std::string osSrcNoData{};

};

static std::unique_ptr<GDALArgumentParser> GDALFootprintAppOptionsGetParser(

    GDALFootprintOptions *psOptions,

    GDALFootprintOptionsForBinary *psOptionsForBinary)

{

    auto argParser = std::make_unique<GDALArgumentParser>(

        "gdal_footprint", /* bForBinary=*/psOptionsForBinary != nullptr);

    argParser->add_description(_("Compute footprint of a raster."));

    argParser->add_epilog(_("For more details, consult "

                            "https://gdal.org/programs/gdal_footprint.html"));

    argParser->add_argument("-b")

        .metavar("<band>")

        .scan<'i', int>()

        .append()

        .store_into(psOptions->anBands)

        .help(_("Band(s) of interest."));

    argParser->add_argument("-combine_bands")

        .choices("union", "intersection")

        .action([psOptions](const auto &s)

                { psOptions->bCombineBandsUnion = s == "union"; })

        .default_value("union")

        .help(_("Defines how the mask bands of the selected bands are combined "

                "to generate a single mask band, before being vectorized."));

    {

        auto &group = argParser->add_mutually_exclusive_group();

        group.add_argument("-ovr")

            .metavar("<index>")

            .scan<'i', int>()

            .store_into(psOptions->nOvrIndex)

            .help(_("Defines which overview level of source file must be used, "

                    "when overviews are available on the source raster."));

        group.add_argument("-srcnodata")

            .metavar("\"<value>[ <value>]...\"")

            .store_into(psOptions->osSrcNoData)

            .help(_("Set nodata value(s) for input bands."));

    }

    argParser->add_argument("-t_cs")

        .choices("pixel", "georef")

        .default_value("georef")

        .action(

            [psOptions](const auto &s)

            {

                const bool georefSet{s == "georef"};

                psOptions->bOutCSGeoref = georefSet;

                psOptions->bOutCSGeorefRequested = georefSet;

            })

        .help(_("Target coordinate system."));

    // Note: no store_into (requires post validation)

    argParser->add_argument("-t_srs")

        .metavar("<srs_def>")

        .help(_("Target CRS of the output file.."));

    argParser->add_argument("-split_polys")

        .flag()

        .store_into(psOptions->bSplitPolys)

        .help(_("Split multipolygons into several features each one with a "

                "single polygon."));

    argParser->add_argument("-convex_hull")

        .flag()

        .store_into(psOptions->bConvexHull)

        .help(_("Compute the convex hull of the (multi)polygons."));

    argParser->add_argument("-densify")

        .metavar("<value>")

        .scan<'g', double>()

        .store_into(psOptions->dfDensifyDistance)

        .help(_("The specified value of this option is the maximum distance "

                "between 2 consecutive points of the output geometry. "));

    argParser->add_argument("-simplify")

        .metavar("<value>")

        .scan<'g', double>()

        .store_into(psOptions->dfSimplifyTolerance)

        .help(_("The specified value of this option is the tolerance used to "

                "merge consecutive points of the output geometry."));

    argParser->add_argument("-min_ring_area")

        .metavar("<value>")

        .scan<'g', double>()

        .store_into(psOptions->dfMinRingArea)

        .help(_("The specified value of this option is the minimum area of a "

                "ring to be considered."));

    // Note: no store_into (requires post validation)

    argParser->add_argument("-max_points")

        .metavar("<value>|unlimited")

        .default_value("100")

        .help(_("The maximum number of points in the output geometry."));

    if (psOptionsForBinary)

    {

        argParser->add_quiet_argument(&psOptionsForBinary->bQuiet);

        argParser->add_open_options_argument(

            psOptionsForBinary->aosOpenOptions);

    }

    argParser->add_output_format_argument(psOptions->osFormat);

    {

        auto &group = argParser->add_mutually_exclusive_group();

        group.add_argument("-location_field_name")

            .metavar("<field_name>")

            .default_value("location")

            .store_into(psOptions->osLocationFieldName)

            .help(_(

                "Specifies the name of the field in the resulting vector "

                "dataset where the path of the input dataset will be stored."));

        group.add_argument("-no_location")

            .flag()

            .store_into(psOptions->bClearLocation)

            .help(

                _("Turns off the writing of the path of the input dataset as a "

                  "field in the output vector dataset."));

    }

    argParser->add_argument("-write_absolute_path")

        .flag()

        .store_into(psOptions->bAbsolutePath)

        .help(_("Enables writing the absolute path of the input dataset."));

    argParser->add_layer_creation_options_argument(psOptions->aosLCO);

    argParser->add_dataset_creation_options_argument(psOptions->aosDSCO);

    argParser->add_argument("-lyr_name")

        .metavar("<value>")

        .store_into(psOptions->osDestLayerName)

        .help(_("Name of the target layer."));

    if (psOptionsForBinary)

    {

        argParser->add_argument("-overwrite")

            .flag()

            .store_into(psOptionsForBinary->bOverwrite)

            .help(_("Overwrite the target layer if it exists."));

        argParser->add_argument("src_filename")

            .metavar("<src_filename>")

            .store_into(psOptionsForBinary->osSource)

            .help(_("Source raster file name."));

        argParser->add_argument("dst_filename")

            .metavar("<dst_filename>")

            .store_into(psOptionsForBinary->osDest)

            .help(_("Destination vector file name."));

    }

    return argParser;

}

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

/*                        GDALFootprintMaskBand                         */

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

class GDALFootprintMaskBand final : public GDALRasterBand

{

    GDALRasterBand *m_poSrcBand = nullptr;

    CPL_DISALLOW_COPY_ASSIGN(GDALFootprintMaskBand)

  public:

    explicit GDALFootprintMaskBand(GDALRasterBand *poSrcBand)

        : m_poSrcBand(poSrcBand)

    {

        nRasterXSize = m_poSrcBand->GetXSize();

        nRasterYSize = m_poSrcBand->GetYSize();

        eDataType = GDT_UInt8;

        m_poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);

    }

  protected:

    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;

    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,

                     int nYSize, void *pData, int nBufXSize, int nBufYSize,

                     GDALDataType eBufType, GSpacing nPixelSpace,

                     GSpacing nLineSpace,

                     GDALRasterIOExtraArg *psExtraArg) override;

};

CPLErr GDALFootprintMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,

                                         void *pData)

{

    int nWindowXSize;

    int nWindowYSize;

    m_poSrcBand->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,

                                    &nWindowYSize);

    GDALRasterIOExtraArg sExtraArg;

    INIT_RASTERIO_EXTRA_ARG(sExtraArg);

    return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,

                     nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,

                     pData, nWindowXSize, nWindowYSize, GDT_UInt8, 1,

                     nBlockXSize, &sExtraArg);

}

CPLErr GDALFootprintMaskBand::IRasterIO(

    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,

    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,

    GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)

{

    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&

        eBufType == GDT_UInt8 && nPixelSpace == 1)

    {

        // Request when band seen as the mask band for GDALPolygonize()

        if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,

                                  nBufXSize, nBufYSize, eBufType, nPixelSpace,

                                  nLineSpace, psExtraArg) != CE_None)

        {

            return CE_Failure;

        }

        GByte *pabyData = static_cast<GByte *>(pData);

        for (int iY = 0; iY < nYSize; ++iY)

        {

            for (int iX = 0; iX < nXSize; ++iX)

            {

                if (pabyData[iX])

                    pabyData[iX] = 1;

            }

            pabyData += nLineSpace;

        }

        return CE_None;

    }

    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&

        eBufType == GDT_Int64 &&

        nPixelSpace == static_cast<int>(sizeof(int64_t)) &&

        (nLineSpace % nPixelSpace) == 0)

    {

        // Request when band seen as the value band for GDALPolygonize()

        if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,

                                  nBufXSize, nBufYSize, eBufType, nPixelSpace,

                                  nLineSpace, psExtraArg) != CE_None)

        {

            return CE_Failure;

        }

        int64_t *panData = static_cast<int64_t *>(pData);

        for (int iY = 0; iY < nYSize; ++iY)

        {

            for (int iX = 0; iX < nXSize; ++iX)

            {

                if (panData[iX])

                    panData[iX] = 1;

            }

            panData += (nLineSpace / nPixelSpace);

        }

        return CE_None;

    }

    return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,

                                     pData, nBufXSize, nBufYSize, eBufType,

                                     nPixelSpace, nLineSpace, psExtraArg);

}

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

/*                    GDALFootprintCombinedMaskBand                     */

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

class GDALFootprintCombinedMaskBand final : public GDALRasterBand

{

    std::vector<GDALRasterBand *> m_apoSrcBands{};

    /** Whether to combine bands unioning (true) or intersecting (false) */

    bool m_bUnion = false;

  public:

    explicit GDALFootprintCombinedMaskBand(

        const std::vector<GDALRasterBand *> &apoSrcBands, bool bUnion)

        : m_apoSrcBands(apoSrcBands), m_bUnion(bUnion)

    {

        nRasterXSize = m_apoSrcBands[0]->GetXSize();

        nRasterYSize = m_apoSrcBands[0]->GetYSize();

        eDataType = GDT_UInt8;

        m_apoSrcBands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);

    }

  protected:

    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;

    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,

                     int nYSize, void *pData, int nBufXSize, int nBufYSize,

                     GDALDataType eBufType, GSpacing nPixelSpace,

                     GSpacing nLineSpace,

                     GDALRasterIOExtraArg *psExtraArg) override;

};

CPLErr GDALFootprintCombinedMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,

                                                 void *pData)

{

    int nWindowXSize;

    int nWindowYSize;

    m_apoSrcBands[0]->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,

                                         &nWindowYSize);

    GDALRasterIOExtraArg sExtraArg;

    INIT_RASTERIO_EXTRA_ARG(sExtraArg);

    return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,

                     nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,

                     pData, nWindowXSize, nWindowYSize, GDT_UInt8, 1,

                     nBlockXSize, &sExtraArg);

}

CPLErr GDALFootprintCombinedMaskBand::IRasterIO(

    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,

    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,

    GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)

{

    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&

        eBufType == GDT_UInt8 && nPixelSpace == 1)

    {

        // Request when band seen as the mask band for GDALPolygonize()

        {

            GByte *pabyData = static_cast<GByte *>(pData);

            for (int iY = 0; iY < nYSize; ++iY)

            {

                memset(pabyData, m_bUnion ? 0 : 1, nXSize);

                pabyData += nLineSpace;

            }

        }

        std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);

        for (auto poBand : m_apoSrcBands)

        {

            if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,

                                 abyTmp.data(), nBufXSize, nBufYSize, GDT_UInt8,

                                 1, nXSize, psExtraArg) != CE_None)

            {

                return CE_Failure;

            }

            GByte *pabyData = static_cast<GByte *>(pData);

            size_t iTmp = 0;

            for (int iY = 0; iY < nYSize; ++iY)

            {

                if (m_bUnion)

                {

                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)

                    {

                        if (abyTmp[iTmp])

                            pabyData[iX] = 1;

                    }

                }

                else

                {

                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)

                    {

                        if (abyTmp[iTmp] == 0)

                            pabyData[iX] = 0;

                    }

                }

                pabyData += nLineSpace;

            }

        }

        return CE_None;

    }

    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&

        eBufType == GDT_Int64 &&

        nPixelSpace == static_cast<int>(sizeof(int64_t)) &&

        (nLineSpace % nPixelSpace) == 0)

    {

        // Request when band seen as the value band for GDALPolygonize()

        {

            int64_t *panData = static_cast<int64_t *>(pData);

            for (int iY = 0; iY < nYSize; ++iY)

            {

                if (m_bUnion)

                {

                    memset(panData, 0, nXSize * sizeof(int64_t));

                }

                else

                {

                    int64_t nOne = 1;

                    GDALCopyWords(&nOne, GDT_Int64, 0, panData, GDT_Int64,

                                  sizeof(int64_t), nXSize);

                }

                panData += (nLineSpace / nPixelSpace);

            }

        }

        std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);

        for (auto poBand : m_apoSrcBands)

        {

            if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,

                                 abyTmp.data(), nBufXSize, nBufYSize, GDT_UInt8,

                                 1, nXSize, psExtraArg) != CE_None)

            {

                return CE_Failure;

            }

            size_t iTmp = 0;

            int64_t *panData = static_cast<int64_t *>(pData);

            for (int iY = 0; iY < nYSize; ++iY)

            {

                if (m_bUnion)

                {

                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)

                    {

                        if (abyTmp[iTmp])

                            panData[iX] = 1;

                    }

                }

                else

                {

                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)

                    {

                        if (abyTmp[iTmp] == 0)

                            panData[iX] = 0;

                    }

                }

                panData += (nLineSpace / nPixelSpace);

            }

        }

        return CE_None;

    }

    return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,

                                     pData, nBufXSize, nBufYSize, eBufType,

                                     nPixelSpace, nLineSpace, psExtraArg);

}

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

/*                    GetOutputLayerAndUpdateDstDS()                    */

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

static OGRLayer *

GetOutputLayerAndUpdateDstDS(const char *pszDest, GDALDatasetH &hDstDS,

                             GDALDataset *poSrcDS,

                             const GDALFootprintOptions *psOptions)

{

    if (pszDest == nullptr)

        pszDest = GDALGetDescription(hDstDS);

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

    /*      Create output dataset if needed                                 */

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

    const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;

    GDALDriverH hDriver = nullptr;

    if (bCreateOutput)

    {

        std::string osFormat(psOptions->osFormat);

        if (osFormat.empty())

        {

            const auto aoDrivers = GetOutputDriversFor(pszDest, GDAL_OF_VECTOR);

            if (aoDrivers.empty())

            {

                CPLError(CE_Failure, CPLE_AppDefined,

                         "Cannot guess driver for %s", pszDest);

                return nullptr;

            }

            else

            {

                if (aoDrivers.size() > 1)

                {

                    CPLError(CE_Warning, CPLE_AppDefined,

                             "Several drivers matching %s extension. Using %s",

                             CPLGetExtensionSafe(pszDest).c_str(),

                             aoDrivers[0].c_str());

                }

                osFormat = aoDrivers[0];

            }

        }

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

        /*      Find the output driver. */

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

        hDriver = GDALGetDriverByName(osFormat.c_str());

        CSLConstList papszDriverMD =

            hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;

        if (hDriver == nullptr ||

            !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_VECTOR,

                                              "FALSE")) ||

            !CPLTestBool(

                CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))

        {

            CPLError(CE_Failure, CPLE_NotSupported,

                     "Output driver `%s' not recognised or does not support "

                     "direct output file creation.",

                     osFormat.c_str());

            return nullptr;

        }

        hDstDS = GDALCreate(hDriver, pszDest, 0, 0, 0, GDT_Unknown,

                            psOptions->aosDSCO.List());

        if (!hDstDS)

        {

            return nullptr;

        }

    }

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

    /*      Open or create target layer.                                    */

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

    auto poDstDS = GDALDataset::FromHandle(hDstDS);

    OGRLayer *poLayer = nullptr;

    if (!bCreateOutput)

    {

        if (poDstDS->GetLayerCount() == 1 && poDstDS->GetDriver() &&

            EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))

        {

            poLayer = poDstDS->GetLayer(0);

        }

        else if (!psOptions->osDestLayerName.empty())

        {

            poLayer =

                poDstDS->GetLayerByName(psOptions->osDestLayerName.c_str());

            if (!poLayer)

            {

                CPLError(CE_Failure, CPLE_AppDefined, "Cannot find layer %s",

                         psOptions->osDestLayerName.c_str());

                return nullptr;

            }

        }

        else

        {

            poLayer = poDstDS->GetLayerByName(DEFAULT_LAYER_NAME);

        }

    }

    if (!poLayer)

    {

        std::string osDestLayerName = psOptions->osDestLayerName;

        if (osDestLayerName.empty())

        {

            if (poDstDS->GetDriver() &&

                EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))

            {

                osDestLayerName = CPLGetBasenameSafe(pszDest);

            }

            else

            {

                osDestLayerName = DEFAULT_LAYER_NAME;

            }

        }

        std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS;

        if (psOptions->bOutCSGeoref)

        {

            if (!psOptions->oOutputSRS.IsEmpty())

            {

                poSRS.reset(psOptions->oOutputSRS.Clone());

            }

            else if (auto poSrcSRS = poSrcDS->GetSpatialRef())

            {

                poSRS.reset(poSrcSRS->Clone());

            }

        }

        poLayer = poDstDS->CreateLayer(

            osDestLayerName.c_str(), poSRS.get(),

            psOptions->bSplitPolys ? wkbPolygon : wkbMultiPolygon,

            const_cast<char **>(psOptions->aosLCO.List()));

        if (!psOptions->osLocationFieldName.empty())

        {

            OGRFieldDefn oFieldDefn(psOptions->osLocationFieldName.c_str(),

                                    OFTString);

            if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)

                return nullptr;

        }

    }

    return poLayer;

}

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

/*                 GeoTransformCoordinateTransformation                 */

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

class GeoTransformCoordinateTransformation final

    : public OGRCoordinateTransformation

{

    const GDALGeoTransform &m_gt;

  public:

    explicit GeoTransformCoordinateTransformation(const GDALGeoTransform &gt)

        : m_gt(gt)

    {

    }

    const OGRSpatialReference *GetSourceCS() const override;

    const OGRSpatialReference *GetTargetCS() const override

    {

        return nullptr;

    }

    OGRCoordinateTransformation *Clone() const override

    {

        return new GeoTransformCoordinateTransformation(m_gt);

    }

    OGRCoordinateTransformation *GetInverse() const override

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "GeoTransformCoordinateTransformation::GetInverse() not "

                 "implemented");

        return nullptr;

    }

    int Transform(size_t nCount, double *x, double *y, double * /* z */,

                  double * /* t */, int *pabSuccess) override

    {

        for (size_t i = 0; i < nCount; ++i)

        {

            const double X = m_gt.xorig + x[i] * m_gt.xscale + y[i] * m_gt.xrot;

            const double Y = m_gt.yorig + x[i] * m_gt.yrot + y[i] * m_gt.yscale;

            x[i] = X;

            y[i] = Y;

            if (pabSuccess)

                pabSuccess[i] = TRUE;

        }

        return TRUE;

    }

};

const OGRSpatialReference *

GeoTransformCoordinateTransformation::GetSourceCS() const

{

    return nullptr;

}

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

/*                            CountPoints()                             */

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

static size_t CountPoints(const OGRGeometry *poGeom)

{

    if (poGeom->getGeometryType() == wkbMultiPolygon)

    {

        size_t n = 0;

        for (auto *poPoly : poGeom->toMultiPolygon())

        {

            n += CountPoints(poPoly);

        }

        return n;

    }

    else if (poGeom->getGeometryType() == wkbPolygon)

    {

        size_t n = 0;

        for (auto *poRing : poGeom->toPolygon())

        {

            n += poRing->getNumPoints() - 1;

        }

        return n;

    }

    return 0;

}

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

/*                   GetMinDistanceBetweenTwoPoints()                   */

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

static double GetMinDistanceBetweenTwoPoints(const OGRGeometry *poGeom)

{

    if (poGeom->getGeometryType() == wkbMultiPolygon)

    {

        double v = std::numeric_limits<double>::max();

        for (auto *poPoly : poGeom->toMultiPolygon())

        {

            v = std::min(v, GetMinDistanceBetweenTwoPoints(poPoly));

        }

        return v;

    }

    else if (poGeom->getGeometryType() == wkbPolygon)

    {

        double v = std::numeric_limits<double>::max();

        for (auto *poRing : poGeom->toPolygon())

        {

            v = std::min(v, GetMinDistanceBetweenTwoPoints(poRing));

        }

        return v;

    }

    else if (poGeom->getGeometryType() == wkbLineString)

    {

        double v = std::numeric_limits<double>::max();

        const auto poLS = poGeom->toLineString();

        const int nNumPoints = poLS->getNumPoints();

        for (int i = 0; i < nNumPoints - 1; ++i)

        {

            const double x1 = poLS->getX(i);

            const double y1 = poLS->getY(i);

            const double x2 = poLS->getX(i + 1);

            const double y2 = poLS->getY(i + 1);

            const double d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);

            if (d > 0)

                v = std::min(v, d);

        }

        return sqrt(v);

    }

    return 0;

}

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

/*                        GDALFootprintProcess()                        */

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

static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,

                                 const GDALFootprintOptions *psOptions)

{

    std::unique_ptr<OGRCoordinateTransformation> poCT_SRS;

    const OGRSpatialReference *poDstSRS = poDstLayer->GetSpatialRef();

    if (!psOptions->oOutputSRS.IsEmpty())

        poDstSRS = &(psOptions->oOutputSRS);

    if (poDstSRS)

    {

        auto poSrcSRS = poSrcDS->GetSpatialRef();

        if (!poSrcSRS)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Output layer has CRS, but input is not georeferenced");

            return false;

        }

        poCT_SRS.reset(OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS));

        if (!poCT_SRS)

            return false;

    }

    std::vector<int> anBands = psOptions->anBands;

    const int nBandCount = poSrcDS->GetRasterCount();

    if (anBands.empty())

    {

        for (int i = 1; i <= nBandCount; ++i)

            anBands.push_back(i);

    }

    std::vector<GDALRasterBand *> apoSrcMaskBands;

    const CPLStringList aosSrcNoData(

        CSLTokenizeString2(psOptions->osSrcNoData.c_str(), " ", 0));

    std::vector<double> adfSrcNoData;

    if (!psOptions->osSrcNoData.empty())

    {

        if (aosSrcNoData.size() != 1 &&

            static_cast<size_t>(aosSrcNoData.size()) != anBands.size())

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Number of values in -srcnodata should be 1 or the number "

                     "of bands");

            return false;

        }

        for (int i = 0; i < aosSrcNoData.size(); ++i)

        {

            adfSrcNoData.emplace_back(CPLAtof(aosSrcNoData[i]));

        }

    }

    bool bGlobalMask = true;

    std::vector<std::unique_ptr<GDALRasterBand>> apoTmpNoDataMaskBands;

    for (size_t i = 0; i < anBands.size(); ++i)

    {

        const int nBand = anBands[i];

        if (nBand <= 0 || nBand > nBandCount)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Invalid band number: %d",

                     nBand);

            return false;

        }

        auto poBand = poSrcDS->GetRasterBand(nBand);

        if (!adfSrcNoData.empty())

        {

            bGlobalMask = false;

            apoTmpNoDataMaskBands.emplace_back(

                std::make_unique<GDALNoDataMaskBand>(

                    poBand, adfSrcNoData.size() == 1 ? adfSrcNoData[0]

                                                     : adfSrcNoData[i]));

            apoSrcMaskBands.push_back(apoTmpNoDataMaskBands.back().get());

        }

        else

        {

            GDALRasterBand *poMaskBand;

            const int nMaskFlags = poBand->GetMaskFlags();

            if (poBand->GetColorInterpretation() == GCI_AlphaBand)

            {