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

 *

 * Project:  GDAL

 * Purpose:  gdal "raster tile" subcommand

 * Author:   Even Rouault <even dot rouault at spatialys.com>

 *

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

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

 *

 * SPDX-License-Identifier: MIT

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

#include "gdalalg_raster_tile.h"

#include "cpl_conv.h"

#include "cpl_json.h"

#include "cpl_mem_cache.h"

#include "cpl_spawn.h"

#include "cpl_time.h"

#include "cpl_vsi_virtual.h"

#include "cpl_worker_thread_pool.h"

#include "gdal_alg_priv.h"

#include "gdal_priv.h"

#include "gdalgetgdalpath.h"

#include "gdalwarper.h"

#include "gdal_utils.h"

#include "ogr_spatialref.h"

#include "memdataset.h"

#include "tilematrixset.hpp"

#include "ogr_p.h"

#include <algorithm>

#include <array>

#include <atomic>

#include <cinttypes>

#include <cmath>

#include <mutex>

#include <utility>

#include <thread>

#ifdef USE_NEON_OPTIMIZATIONS

#include "include_sse2neon.h"

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

#include <emmintrin.h>

#if defined(__SSSE3__) || defined(__AVX__)

#include <tmmintrin.h>

#endif

#if defined(__SSE4_1__) || defined(__AVX__)

#include <smmintrin.h>

#endif

#endif

#if defined(__x86_64) || defined(_M_X64) || defined(USE_NEON_OPTIMIZATIONS)

#define USE_PAETH_SSE2

#endif

#ifndef _WIN32

#define FORK_ALLOWED

#endif

#include "cpl_zlib_header.h"  // for crc32()

//! @cond Doxygen_Suppress

#ifndef _

#define _(x) (x)

#endif

// Unlikely substring to appear in stdout. We do that in case some GDAL

// driver would output on stdout.

constexpr const char PROGRESS_MARKER[] = {'!', '.', 'x'};

constexpr const char END_MARKER[] = {'?', 'E', '?', 'N', '?', 'D', '?'};

constexpr const char ERROR_START_MARKER[] = {'%', 'E', '%', 'R', '%', 'R',

                                             '%', '_', '%', 'S', '%', 'T',

                                             '%', 'A', '%', 'R', '%', 'T'};

constexpr const char *STOP_MARKER = "STOP\n";

namespace

{

struct BandMetadata

{

    std::string osDescription{};

    GDALDataType eDT{};

    GDALColorInterp eColorInterp{};

    std::string osCenterWaveLength{};

    std::string osFWHM{};

};

}  // namespace

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

/*                     GetThresholdMinTilesPerJob()                     */

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

static int GetThresholdMinThreadsForSpawn()

{

    // Minimum number of threads for automatic switch to spawning

    constexpr int THRESHOLD_MIN_THREADS_FOR_SPAWN = 8;

    // Config option for test only

    return std::max(1, atoi(CPLGetConfigOption(

                           "GDAL_THRESHOLD_MIN_THREADS_FOR_SPAWN",

                           CPLSPrintf("%d", THRESHOLD_MIN_THREADS_FOR_SPAWN))));

}

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

/*                     GetThresholdMinTilesPerJob()                     */

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

static int GetThresholdMinTilesPerJob()

{

    // Minimum number of tiles per job to decide for automatic switch to spawning

    constexpr int THRESHOLD_TILES_PER_JOB = 100;

    // Config option for test only

    return std::max(

        1, atoi(CPLGetConfigOption("GDAL_THRESHOLD_MIN_TILES_PER_JOB",

                                   CPLSPrintf("%d", THRESHOLD_TILES_PER_JOB))));

}

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

/*          GDALRasterTileAlgorithm::GDALRasterTileAlgorithm()          */

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

GDALRasterTileAlgorithm::GDALRasterTileAlgorithm(bool standaloneStep)

    : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,

                                      ConstructorOptions()

                                          .SetStandaloneStep(standaloneStep)

                                          .SetInputDatasetMaxCount(1)

                                          .SetAddDefaultArguments(false)

                                          .SetInputDatasetAlias("dataset"))

{

    if (standaloneStep)

        AddProgressArg();

    AddArg("spawned", 0, _("Whether this is a spawned worker"),

           &m_spawned)

        .SetHidden();  // Used in spawn mode

#ifdef FORK_ALLOWED

    AddArg("forked", 0, _("Whether this is a forked worker"),

           &m_forked)

        .SetHidden();  // Used in forked mode

#else

    CPL_IGNORE_RET_VAL(m_forked);

#endif

    AddArg("config-options-in-stdin", 0, _(""), &m_dummy)

        .SetHidden();  // Used in spawn mode

    AddArg("ovr-zoom-level", 0, _("Overview zoom level to compute"),

           &m_ovrZoomLevel)

        .SetMinValueIncluded(0)

        .SetHidden();  // Used in spawn mode

    AddArg("ovr-min-x", 0, _("Minimum tile X coordinate"), &m_minOvrTileX)

        .SetMinValueIncluded(0)

        .SetHidden();  // Used in spawn mode

    AddArg("ovr-max-x", 0, _("Maximum tile X coordinate"), &m_maxOvrTileX)

        .SetMinValueIncluded(0)

        .SetHidden();  // Used in spawn mode

    AddArg("ovr-min-y", 0, _("Minimum tile Y coordinate"), &m_minOvrTileY)

        .SetMinValueIncluded(0)

        .SetHidden();  // Used in spawn mode

    AddArg("ovr-max-y", 0, _("Maximum tile Y coordinate"), &m_maxOvrTileY)

        .SetMinValueIncluded(0)

        .SetHidden();  // Used in spawn mode

    if (standaloneStep)

    {

        AddRasterInputArgs(/* openForMixedRasterVector = */ false,

                           /* hiddenForCLI = */ false);

    }

    else

    {

        AddRasterHiddenInputDatasetArg();

    }

    m_format = "PNG";

    AddOutputFormatArg(&m_format)

        .SetDefault(m_format)

        .AddMetadataItem(

            GAAMDI_REQUIRED_CAPABILITIES,

            {GDAL_DCAP_RASTER, GDAL_DCAP_CREATECOPY, GDAL_DMD_EXTENSIONS})

        .AddMetadataItem(GAAMDI_VRT_COMPATIBLE, {"false"});

    AddCreationOptionsArg(&m_creationOptions);

    AddArg(GDAL_ARG_NAME_OUTPUT, 'o', _("Output directory"), &m_outputDir)

        .SetRequired()

        .SetIsInput()

        .SetMinCharCount(1)

        .SetPositional();

    std::vector<std::string> tilingSchemes{"raster"};

    for (const std::string &scheme :

         gdal::TileMatrixSet::listPredefinedTileMatrixSets(/* hidden = */ true))

    {

        auto poTMS = gdal::TileMatrixSet::parse(scheme.c_str());

        OGRSpatialReference oSRS_TMS;

        if (poTMS && !poTMS->hasVariableMatrixWidth() &&

            oSRS_TMS.SetFromUserInput(poTMS->crs().c_str()) == OGRERR_NONE)

        {

            std::string identifier = scheme == "GoogleMapsCompatible"

                                         ? "WebMercatorQuad"

                                         : poTMS->identifier();

            m_mapTileMatrixIdentifierToScheme[identifier] = scheme;

            tilingSchemes.push_back(std::move(identifier));

        }

    }

    AddArg("tiling-scheme", 0, _("Tiling scheme"), &m_tilingScheme)

        .SetDefault("WebMercatorQuad")

        .SetChoices(tilingSchemes)

        .SetHiddenChoices(

            "GoogleMapsCompatible",  // equivalent of WebMercatorQuad

            "mercator",              // gdal2tiles equivalent of WebMercatorQuad

            "GlobalGeodeticOriginLat270"  // gdal2tiles geodetic without --tmscompatible

        );

    AddArg("min-zoom", 0, _("Minimum zoom level"), &m_minZoomLevel)

        .SetMinValueIncluded(0);

    // Only used by PMTiles driver for now

    AddArg("min-zoom-single-tile", 0,

           _("Determine minimum zoom level to produce a single tile"),

           &m_minZoomLevelSingleTile)

        .SetHidden();

    AddArg("max-zoom", 0, _("Maximum zoom level"), &m_maxZoomLevel)

        .SetMinValueIncluded(0);

    AddArg("min-x", 0, _("Minimum tile X coordinate"), &m_minTileX)

        .SetMinValueIncluded(0);

    AddArg("max-x", 0, _("Maximum tile X coordinate"), &m_maxTileX)

        .SetMinValueIncluded(0);

    AddArg("min-y", 0, _("Minimum tile Y coordinate"), &m_minTileY)

        .SetMinValueIncluded(0);

    AddArg("max-y", 0, _("Maximum tile Y coordinate"), &m_maxTileY)

        .SetMinValueIncluded(0);

    AddArg("no-intersection-ok", 0,

           _("Whether dataset extent not intersecting tile matrix is only a "

             "warning"),

           &m_noIntersectionIsOK);

    AddArg("resampling", 'r', _("Resampling method for max zoom"),

           &m_resampling)

        .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",

                    "average", "rms", "mode", "min", "max", "med", "q1", "q3",

                    "sum")

        .SetDefault("cubic")

        .SetHiddenChoices("near");

    AddArg("overview-resampling", 0, _("Resampling method for overviews"),

           &m_overviewResampling)

        .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",

                    "average", "rms", "mode", "min", "max", "med", "q1", "q3",

                    "sum")

        .SetHiddenChoices("near");

    AddArg("convention", 0,

           _("Tile numbering convention: xyz (from top) or tms (from bottom)"),

           &m_convention)

        .SetDefault(m_convention)

        .SetChoices("xyz", "tms");

    AddArg("tile-size", 0, _("Override default tile size"), &m_tileSize)

        .SetMinValueIncluded(64)

        .SetMaxValueIncluded(32768);

    AddArg("add-alpha", 0, _("Whether to force adding an alpha channel"),

           &m_addalpha)

        .SetMutualExclusionGroup("alpha");

    AddArg("no-alpha", 0, _("Whether to disable adding an alpha channel"),

           &m_noalpha)

        .SetMutualExclusionGroup("alpha");

    auto &dstNoDataArg =

        AddArg("output-nodata", 0, _("Output nodata value"), &m_dstNoData)

            .AddHiddenAlias("dst-nodata");

    AddArg("skip-blank", 0, _("Do not generate blank tiles"), &m_skipBlank);

    {

        auto &arg = AddArg("metadata", 0,

                           _("Add metadata item to output tiles"), &m_metadata)

                        .SetMetaVar("<KEY>=<VALUE>")

                        .SetPackedValuesAllowed(false);

        arg.AddValidationAction([this, &arg]()

                                { return ParseAndValidateKeyValue(arg); });

        arg.AddHiddenAlias("mo");

    }

    AddArg("copy-src-metadata", 0,

           _("Whether to copy metadata from source dataset"),

           &m_copySrcMetadata);

    AddArg("aux-xml", 0, _("Generate .aux.xml sidecar files when needed"),

           &m_auxXML);

    AddArg("kml", 0, _("Generate KML files"), &m_kml);

    AddArg("resume", 0, _("Generate only missing files"), &m_resume);

    AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);

    AddArg("parallel-method", 0,

#ifdef FORK_ALLOWED

           _("Parallelization method (thread, spawn, fork)")

#else

           _("Parallelization method (thread / spawn)")

#endif

               ,

           &m_parallelMethod)

        .SetChoices("thread", "spawn"

#ifdef FORK_ALLOWED

                    ,

                    "fork"

#endif

        );

    constexpr const char *ADVANCED_RESAMPLING_CATEGORY = "Advanced Resampling";

    auto &excludedValuesArg =

        AddArg("excluded-values", 0,

               _("Tuples of values (e.g. <R>,<G>,<B> or (<R1>,<G1>,<B1>),"

                 "(<R2>,<G2>,<B2>)) that must beignored as contributing source "

                 "pixels during (average) resampling"),

               &m_excludedValues)

            .SetCategory(ADVANCED_RESAMPLING_CATEGORY);

    auto &excludedValuesPctThresholdArg =

        AddArg(

            "excluded-values-pct-threshold", 0,

            _("Minimum percentage of source pixels that must be set at one of "

              "the --excluded-values to cause the excluded value to be used as "

              "the target pixel value"),

            &m_excludedValuesPctThreshold)

            .SetDefault(m_excludedValuesPctThreshold)

            .SetMinValueIncluded(0)

            .SetMaxValueIncluded(100)

            .SetCategory(ADVANCED_RESAMPLING_CATEGORY);

    auto &nodataValuesPctThresholdArg =

        AddArg(

            "nodata-values-pct-threshold", 0,

            _("Minimum percentage of source pixels that must be set at one of "

              "nodata (or alpha=0 or any other way to express transparent pixel"

              "to cause the target pixel value to be transparent"),

            &m_nodataValuesPctThreshold)

            .SetDefault(m_nodataValuesPctThreshold)

            .SetMinValueIncluded(0)

            .SetMaxValueIncluded(100)

            .SetCategory(ADVANCED_RESAMPLING_CATEGORY);

    constexpr const char *PUBLICATION_CATEGORY = "Publication";

    AddArg("webviewer", 0, _("Web viewer to generate"), &m_webviewers)

        .SetDefault("all")

        .SetChoices("none", "all", "leaflet", "openlayers", "mapml", "stac")

        .SetCategory(PUBLICATION_CATEGORY);

    AddArg("url", 0,

           _("URL address where the generated tiles are going to be published"),

           &m_url)

        .SetCategory(PUBLICATION_CATEGORY);

    AddArg("title", 0, _("Title of the map"), &m_title)

        .SetCategory(PUBLICATION_CATEGORY);

    AddArg("copyright", 0, _("Copyright for the map"), &m_copyright)

        .SetCategory(PUBLICATION_CATEGORY);

    AddArg("mapml-template", 0,

           _("Filename of a template mapml file where variables will be "

             "substituted"),

           &m_mapmlTemplate)

        .SetMinCharCount(1)

        .SetCategory(PUBLICATION_CATEGORY);

    AddValidationAction(

        [this, &dstNoDataArg, &excludedValuesArg,

         &excludedValuesPctThresholdArg, &nodataValuesPctThresholdArg]()

        {

            if (m_minTileX >= 0 && m_maxTileX >= 0 && m_minTileX > m_maxTileX)

            {

                ReportError(CE_Failure, CPLE_IllegalArg,

                            "'min-x' must be lesser or equal to 'max-x'");

                return false;

            }

            if (m_minTileY >= 0 && m_maxTileY >= 0 && m_minTileY > m_maxTileY)

            {

                ReportError(CE_Failure, CPLE_IllegalArg,

                            "'min-y' must be lesser or equal to 'max-y'");

                return false;

            }

            if (m_minZoomLevel >= 0 && m_maxZoomLevel >= 0 &&

                m_minZoomLevel > m_maxZoomLevel)

            {

                ReportError(CE_Failure, CPLE_IllegalArg,

                            "'min-zoom' must be lesser or equal to 'max-zoom'");

                return false;

            }

            if (m_addalpha && dstNoDataArg.IsExplicitlySet())

            {

                ReportError(

                    CE_Failure, CPLE_IllegalArg,

                    "'add-alpha' and 'output-nodata' are mutually exclusive");

                return false;

            }

            for (const auto *arg :

                 {&excludedValuesArg, &excludedValuesPctThresholdArg,

                  &nodataValuesPctThresholdArg})

            {

                if (arg->IsExplicitlySet() && m_resampling != "average")

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "'%s' can only be specified if 'resampling' is "

                                "set to 'average'",

                                arg->GetName().c_str());

                    return false;

                }

                if (arg->IsExplicitlySet() && !m_overviewResampling.empty() &&

                    m_overviewResampling != "average")

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "'%s' can only be specified if "

                                "'overview-resampling' is set to 'average'",

                                arg->GetName().c_str());

                    return false;

                }

            }

            return true;

        });

}

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

/*                      ~GDALRasterTileAlgorithm()                      */

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

GDALRasterTileAlgorithm::~GDALRasterTileAlgorithm()

{

    if (m_poSrcOvrDS)

    {

        m_poSrcOvrDS->ReleaseRef();

    }

}

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

/*                           GetTileIndices()                           */

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

static bool GetTileIndices(gdal::TileMatrixSet::TileMatrix &tileMatrix,

                           bool bInvertAxisTMS, int tileSize,

                           const double adfExtent[4], int &nMinTileX,

                           int &nMinTileY, int &nMaxTileX, int &nMaxTileY,

                           bool noIntersectionIsOK, bool &bIntersects,

                           bool checkRasterOverflow = true)

{

    if (tileSize > 0)

    {

        tileMatrix.mResX *=

            static_cast<double>(tileMatrix.mTileWidth) / tileSize;

        tileMatrix.mResY *=

            static_cast<double>(tileMatrix.mTileHeight) / tileSize;

        tileMatrix.mTileWidth = tileSize;

        tileMatrix.mTileHeight = tileSize;

    }

    if (bInvertAxisTMS)

        std::swap(tileMatrix.mTopLeftX, tileMatrix.mTopLeftY);

    const double dfTileWidth = tileMatrix.mResX * tileMatrix.mTileWidth;

    const double dfTileHeight = tileMatrix.mResY * tileMatrix.mTileHeight;

    constexpr double EPSILON = 1e-3;

    const double dfMinTileX =

        (adfExtent[0] - tileMatrix.mTopLeftX) / dfTileWidth;

    nMinTileX = static_cast<int>(

        std::clamp(std::floor(dfMinTileX + EPSILON), 0.0,

                   static_cast<double>(tileMatrix.mMatrixWidth - 1)));

    const double dfMinTileY =

        (tileMatrix.mTopLeftY - adfExtent[3]) / dfTileHeight;

    nMinTileY = static_cast<int>(

        std::clamp(std::floor(dfMinTileY + EPSILON), 0.0,

                   static_cast<double>(tileMatrix.mMatrixHeight - 1)));

    const double dfMaxTileX =

        (adfExtent[2] - tileMatrix.mTopLeftX) / dfTileWidth;

    nMaxTileX = static_cast<int>(

        std::clamp(std::floor(dfMaxTileX + EPSILON), 0.0,

                   static_cast<double>(tileMatrix.mMatrixWidth - 1)));

    const double dfMaxTileY =

        (tileMatrix.mTopLeftY - adfExtent[1]) / dfTileHeight;

    nMaxTileY = static_cast<int>(

        std::clamp(std::floor(dfMaxTileY + EPSILON), 0.0,

                   static_cast<double>(tileMatrix.mMatrixHeight - 1)));

    bIntersects = (dfMinTileX <= tileMatrix.mMatrixWidth && dfMaxTileX >= 0 &&

                   dfMinTileY <= tileMatrix.mMatrixHeight && dfMaxTileY >= 0);

    if (!bIntersects)

    {

        CPLDebug("gdal_raster_tile",

                 "dfMinTileX=%g dfMinTileY=%g dfMaxTileX=%g dfMaxTileY=%g",

                 dfMinTileX, dfMinTileY, dfMaxTileX, dfMaxTileY);

        CPLError(noIntersectionIsOK ? CE_Warning : CE_Failure, CPLE_AppDefined,

                 "Extent of source dataset is not compatible with extent of "

                 "tile matrix %s",

                 tileMatrix.mId.c_str());

        return noIntersectionIsOK;

    }

    if (checkRasterOverflow)

    {

        if (nMaxTileX - nMinTileX + 1 > INT_MAX / tileMatrix.mTileWidth ||

            nMaxTileY - nMinTileY + 1 > INT_MAX / tileMatrix.mTileHeight)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Too large zoom level");

            return false;

        }

    }

    return true;

}

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

/*                              GetFileY()                              */

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

static int GetFileY(int iY, const gdal::TileMatrixSet::TileMatrix &tileMatrix,

                    const std::string &convention)

{

    return convention == "xyz" ? iY : tileMatrix.mMatrixHeight - 1 - iY;

}

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

/*                            GenerateTile()                            */

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

// Cf http://www.libpng.org/pub/png/spec/1.2/PNG-Filters.html

// for specification of SUB and AVG filters

inline GByte PNG_SUB(int nVal, int nValPrev)

{

    return static_cast<GByte>((nVal - nValPrev) & 0xff);

}

inline GByte PNG_AVG(int nVal, int nValPrev, int nValUp)

{

    return static_cast<GByte>((nVal - (nValPrev + nValUp) / 2) & 0xff);

}

inline GByte PNG_PAETH(int nVal, int nValPrev, int nValUp, int nValUpPrev)

{

    const int p = nValPrev + nValUp - nValUpPrev;

    const int pa = std::abs(p - nValPrev);

    const int pb = std::abs(p - nValUp);

    const int pc = std::abs(p - nValUpPrev);

    if (pa <= pb && pa <= pc)

        return static_cast<GByte>((nVal - nValPrev) & 0xff);

    else if (pb <= pc)

        return static_cast<GByte>((nVal - nValUp) & 0xff);

    else

        return static_cast<GByte>((nVal - nValUpPrev) & 0xff);

}

#ifdef USE_PAETH_SSE2

static inline __m128i abs_epi16(__m128i x)

{

#if defined(__SSSE3__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)

    return _mm_abs_epi16(x);

#else

    __m128i mask = _mm_srai_epi16(x, 15);

    return _mm_sub_epi16(_mm_xor_si128(x, mask), mask);

#endif

}

static inline __m128i blendv(__m128i a, __m128i b, __m128i mask)

{

#if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)

    return _mm_blendv_epi8(a, b, mask);

#else

    return _mm_or_si128(_mm_andnot_si128(mask, a), _mm_and_si128(mask, b));

#endif

}

static inline __m128i PNG_PAETH_SSE2(__m128i up_prev, __m128i up, __m128i prev,

                                     __m128i cur, __m128i &cost)

{

    auto cur_lo = _mm_unpacklo_epi8(cur, _mm_setzero_si128());

    auto prev_lo = _mm_unpacklo_epi8(prev, _mm_setzero_si128());

    auto up_lo = _mm_unpacklo_epi8(up, _mm_setzero_si128());

    auto up_prev_lo = _mm_unpacklo_epi8(up_prev, _mm_setzero_si128());

    auto cur_hi = _mm_unpackhi_epi8(cur, _mm_setzero_si128());

    auto prev_hi = _mm_unpackhi_epi8(prev, _mm_setzero_si128());

    auto up_hi = _mm_unpackhi_epi8(up, _mm_setzero_si128());

    auto up_prev_hi = _mm_unpackhi_epi8(up_prev, _mm_setzero_si128());

    auto pa_lo = _mm_sub_epi16(up_lo, up_prev_lo);

    auto pb_lo = _mm_sub_epi16(prev_lo, up_prev_lo);

    auto pc_lo = _mm_add_epi16(pa_lo, pb_lo);

    pa_lo = abs_epi16(pa_lo);

    pb_lo = abs_epi16(pb_lo);

    pc_lo = abs_epi16(pc_lo);

    auto min_lo = _mm_min_epi16(_mm_min_epi16(pa_lo, pb_lo), pc_lo);

    auto res_lo = blendv(up_prev_lo, up_lo, _mm_cmpeq_epi16(min_lo, pb_lo));

    res_lo = blendv(res_lo, prev_lo, _mm_cmpeq_epi16(min_lo, pa_lo));

    res_lo = _mm_and_si128(_mm_sub_epi16(cur_lo, res_lo), _mm_set1_epi16(0xFF));

    auto cost_lo = blendv(_mm_sub_epi16(_mm_set1_epi16(256), res_lo), res_lo,

                          _mm_cmplt_epi16(res_lo, _mm_set1_epi16(128)));

    auto pa_hi = _mm_sub_epi16(up_hi, up_prev_hi);

    auto pb_hi = _mm_sub_epi16(prev_hi, up_prev_hi);

    auto pc_hi = _mm_add_epi16(pa_hi, pb_hi);

    pa_hi = abs_epi16(pa_hi);

    pb_hi = abs_epi16(pb_hi);

    pc_hi = abs_epi16(pc_hi);

    auto min_hi = _mm_min_epi16(_mm_min_epi16(pa_hi, pb_hi), pc_hi);

    auto res_hi = blendv(up_prev_hi, up_hi, _mm_cmpeq_epi16(min_hi, pb_hi));

    res_hi = blendv(res_hi, prev_hi, _mm_cmpeq_epi16(min_hi, pa_hi));

    res_hi = _mm_and_si128(_mm_sub_epi16(cur_hi, res_hi), _mm_set1_epi16(0xFF));

    auto cost_hi = blendv(_mm_sub_epi16(_mm_set1_epi16(256), res_hi), res_hi,

                          _mm_cmplt_epi16(res_hi, _mm_set1_epi16(128)));

    cost_lo = _mm_add_epi16(cost_lo, cost_hi);

    cost =

        _mm_add_epi32(cost, _mm_unpacklo_epi16(cost_lo, _mm_setzero_si128()));

    cost =

        _mm_add_epi32(cost, _mm_unpackhi_epi16(cost_lo, _mm_setzero_si128()));

    return _mm_packus_epi16(res_lo, res_hi);

}

static int RunPaeth(const GByte *srcBuffer, int nBands,

                    int nSrcBufferBandStride, GByte *outBuffer, int W,

                    int &costPaeth)

{

    __m128i xmm_cost = _mm_setzero_si128();

    int i = 1;

    for (int k = 0; k < nBands; ++k)

    {

        for (i = 1; i + 15 < W; i += 16)

        {

            auto up_prev = _mm_loadu_si128(

                reinterpret_cast<const __m128i *>(srcBuffer - W + (i - 1)));

            auto up = _mm_loadu_si128(

                reinterpret_cast<const __m128i *>(srcBuffer - W + i));

            auto prev = _mm_loadu_si128(

                reinterpret_cast<const __m128i *>(srcBuffer + (i - 1)));

            auto cur = _mm_loadu_si128(

                reinterpret_cast<const __m128i *>(srcBuffer + i));

            auto res = PNG_PAETH_SSE2(up_prev, up, prev, cur, xmm_cost);

            _mm_storeu_si128(reinterpret_cast<__m128i *>(outBuffer + k * W + i),

                             res);

        }

        srcBuffer += nSrcBufferBandStride;

    }

    int32_t ar_cost[4];

    _mm_storeu_si128(reinterpret_cast<__m128i *>(ar_cost), xmm_cost);

    for (int k = 0; k < 4; ++k)

        costPaeth += ar_cost[k];

    return i;

}

#endif  // USE_PAETH_SSE2

static bool GenerateTile(

    GDALDataset *poSrcDS, GDALDriver *m_poDstDriver, const char *pszExtension,

    CSLConstList creationOptions, GDALWarpOperation &oWO,

    const OGRSpatialReference &oSRS_TMS, GDALDataType eWorkingDataType,

    const gdal::TileMatrixSet::TileMatrix &tileMatrix,

    const std::string &outputDirectory, int nBands, const double *pdfDstNoData,

    int nZoomLevel, int iX, int iY, const std::string &convention,

    int nMinTileX, int nMinTileY, bool bSkipBlank, bool bUserAskedForAlpha,

    bool bAuxXML, bool bResume, const std::vector<std::string> &metadata,

    const GDALColorTable *poColorTable, std::vector<GByte> &dstBuffer,

    std::vector<GByte> &tmpBuffer)

{

    const std::string osDirZ = CPLFormFilenameSafe(

        outputDirectory.c_str(), CPLSPrintf("%d", nZoomLevel), nullptr);

    const std::string osDirX =

        CPLFormFilenameSafe(osDirZ.c_str(), CPLSPrintf("%d", iX), nullptr);

    const int iFileY = GetFileY(iY, tileMatrix, convention);

    const std::string osFilename = CPLFormFilenameSafe(

        osDirX.c_str(), CPLSPrintf("%d", iFileY), pszExtension);

    if (bResume)

    {

        VSIStatBufL sStat;

        if (VSIStatL(osFilename.c_str(), &sStat) == 0)

            return true;

    }

    const int nDstXOff = (iX - nMinTileX) * tileMatrix.mTileWidth;

    const int nDstYOff = (iY - nMinTileY) * tileMatrix.mTileHeight;

    memset(dstBuffer.data(), 0, dstBuffer.size());

    const CPLErr eErr = oWO.WarpRegionToBuffer(

        nDstXOff, nDstYOff, tileMatrix.mTileWidth, tileMatrix.mTileHeight,

        dstBuffer.data(), eWorkingDataType);

    if (eErr != CE_None)

        return false;

    bool bDstHasAlpha =

        nBands > poSrcDS->GetRasterCount() ||

        (nBands == poSrcDS->GetRasterCount() &&

         poSrcDS->GetRasterBand(nBands)->GetColorInterpretation() ==

             GCI_AlphaBand);

    const size_t nBytesPerBand = static_cast<size_t>(tileMatrix.mTileWidth) *

                                 tileMatrix.mTileHeight *

                                 GDALGetDataTypeSizeBytes(eWorkingDataType);

    if (bDstHasAlpha && bSkipBlank)

    {

        bool bBlank = true;

        for (size_t i = 0; i < nBytesPerBand && bBlank; ++i)

        {

            bBlank = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 0);

        }

        if (bBlank)

            return true;

    }

    if (bDstHasAlpha && !bUserAskedForAlpha)

    {

        bool bAllOpaque = true;

        for (size_t i = 0; i < nBytesPerBand && bAllOpaque; ++i)

        {

            bAllOpaque = (dstBuffer[(nBands - 1) * nBytesPerBand + i] == 255);

        }

        if (bAllOpaque)

        {

            bDstHasAlpha = false;

            nBands--;

        }

    }

    VSIMkdir(osDirZ.c_str(), 0755);

    VSIMkdir(osDirX.c_str(), 0755);

    const bool bSupportsCreateOnlyVisibleAtCloseTime =

        m_poDstDriver->GetMetadataItem(

            GDAL_DCAP_CREATE_ONLY_VISIBLE_AT_CLOSE_TIME) != nullptr;

    const std::string osTmpFilename = bSupportsCreateOnlyVisibleAtCloseTime

                                          ? osFilename

                                          : osFilename + ".tmp." + pszExtension;

    const int W = tileMatrix.mTileWidth;

    const int H = tileMatrix.mTileHeight;

    constexpr int EXTRA_BYTE_PER_ROW = 1;  // for filter type

    constexpr int EXTRA_ROWS = 2;          // for paethBuffer and paethBufferTmp

    if (!bAuxXML && EQUAL(pszExtension, "png") &&

        eWorkingDataType == GDT_UInt8 && poColorTable == nullptr &&

        pdfDstNoData == nullptr && W <= INT_MAX / nBands &&

        nBands * W <= INT_MAX - EXTRA_BYTE_PER_ROW &&

        H <= INT_MAX - EXTRA_ROWS &&

        EXTRA_BYTE_PER_ROW + nBands * W <= INT_MAX / (H + EXTRA_ROWS) &&

        CSLCount(creationOptions) == 0 &&

        CPLTestBool(

            CPLGetConfigOption("GDAL_RASTER_TILE_USE_PNG_OPTIM", "YES")))

    {

        // This is an optimized code path completely shortcircuiting libpng

        // We manually generate the PNG file using the Average or PAETH filter

        // and ZLIB compressing the whole buffer, hopefully with libdeflate.

        const int nDstBytesPerRow = EXTRA_BYTE_PER_ROW + nBands * W;

        const int nBPB = static_cast<int>(nBytesPerBand);

        bool bBlank = false;

        if (bDstHasAlpha)

        {

            bBlank = true;

            for (int i = 0; i < nBPB && bBlank; ++i)

            {

                bBlank = (dstBuffer[(nBands - 1) * nBPB + i] == 0);

            }

        }

        constexpr GByte PNG_FILTER_SUB = 1;  // horizontal diff

        constexpr GByte PNG_FILTER_AVG = 3;  // average with pixel before and up

        constexpr GByte PNG_FILTER_PAETH = 4;

        if (bBlank)

            tmpBuffer.clear();

        const int tmpBufferSize = cpl::fits_on<int>(nDstBytesPerRow * H);

        try

        {

            // cppcheck-suppress integerOverflowCond

            tmpBuffer.resize(tmpBufferSize + EXTRA_ROWS * nDstBytesPerRow);

        }

        catch (const std::exception &)

        {

            CPLError(CE_Failure, CPLE_OutOfMemory,

                     "Out of memory allocating temporary buffer");

            return false;

        }

        GByte *const paethBuffer = tmpBuffer.data() + tmpBufferSize;

#ifdef USE_PAETH_SSE2

        GByte *const paethBufferTmp =

            tmpBuffer.data() + tmpBufferSize + nDstBytesPerRow;

#endif

        const char *pszGDAL_RASTER_TILE_PNG_FILTER =

            CPLGetConfigOption("GDAL_RASTER_TILE_PNG_FILTER", "");

        const bool bForcePaeth = EQUAL(pszGDAL_RASTER_TILE_PNG_FILTER, "PAETH");

        const bool bForceAvg = EQUAL(pszGDAL_RASTER_TILE_PNG_FILTER, "AVERAGE");

        for (int j = 0; !bBlank && j < H; ++j)

        {

            if (j > 0)

            {

                tmpBuffer[cpl::fits_on<int>(j * nDstBytesPerRow)] =

                    PNG_FILTER_AVG;

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

                {

                    tmpBuffer[1 + j * nDstBytesPerRow + i] =

                        PNG_AVG(dstBuffer[i * nBPB + j * W], 0,

                                dstBuffer[i * nBPB + (j - 1) * W]);

                }

            }

            else

            {

                tmpBuffer[cpl::fits_on<int>(j * nDstBytesPerRow)] =

                    PNG_FILTER_SUB;

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

                {

                    tmpBuffer[1 + j * nDstBytesPerRow + i] =

                        dstBuffer[i * nBPB + j * W];

                }

            }

            if (nBands == 1)

            {

                if (j > 0)

                {

                    int costAvg = 0;

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

                    {

                        const GByte v =

                            PNG_AVG(dstBuffer[0 * nBPB + j * W + i],

                                    dstBuffer[0 * nBPB + j * W + i - 1],

                                    dstBuffer[0 * nBPB + (j - 1) * W + i]);

                        tmpBuffer[1 + j * nDstBytesPerRow + i * nBands + 0] = v;

                        costAvg += (v < 128) ? v : 256 - v;

                    }

                    if (!bForceAvg)

                    {

                        int costPaeth = 0;

                        {

                            const int i = 0;

                            const GByte v = PNG_PAETH(

                                dstBuffer[0 * nBPB + j * W + i], 0,

                                dstBuffer[0 * nBPB + (j - 1) * W + i], 0);

                            paethBuffer[i] = v;

                            costPaeth += (v < 128) ? v : 256 - v;

                        }

#ifdef USE_PAETH_SSE2

                        const int iLimitSSE2 =

                            RunPaeth(dstBuffer.data() + j * W, nBands, nBPB,

                                     paethBuffer, W, costPaeth);

                        int i = iLimitSSE2;

#else

                        int i = 1;

#endif

                        for (; i < W && (costPaeth < costAvg || bForcePaeth);

                             ++i)

                        {

                            const GByte v = PNG_PAETH(

                                dstBuffer[0 * nBPB + j * W + i],

                                dstBuffer[0 * nBPB + j * W + i - 1],

                                dstBuffer[0 * nBPB + (j - 1) * W + i],

                                dstBuffer[0 * nBPB + (j - 1) * W + i - 1]);

                            paethBuffer[i] = v;

                            costPaeth += (v < 128) ? v : 256 - v;

                        }

                        if (costPaeth < costAvg || bForcePaeth)

                        {

                            GByte *out = tmpBuffer.data() +

                                         cpl::fits_on<int>(j * nDstBytesPerRow);

                            *out = PNG_FILTER_PAETH;

                            ++out;

                            memcpy(out, paethBuffer, nDstBytesPerRow - 1);

                        }

                    }

                }

                else

                {

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

                    {

                        tmpBuffer[1 + j * nDstBytesPerRow + i * nBands + 0] =

                            PNG_SUB(dstBuffer[0 * nBPB + j * W + i],

                                    dstBuffer[0 * nBPB + j * W + i - 1]);

                    }

                }

            }

            else if (nBands == 2)

            {

                if (j > 0)