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

 *

 * Project:  GDAL

 * Purpose:  "blend" step of "raster pipeline"

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

 *

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

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

 * Copyright (c) 2009, Frank Warmerdam

 * SPDX-License-Identifier: MIT

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

#include "gdalalg_raster_blend.h"

#include "cpl_conv.h"

#include "gdal_priv.h"

#include "gdal_utils.h"

#include <algorithm>

#include <array>

#include <cassert>

#include <limits>

#include <utility>

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

#define HAVE_SSE2

#include <immintrin.h>

#endif

#ifdef HAVE_SSE2

#include "gdalsse_priv.h"

#endif

//! @cond Doxygen_Suppress

#ifndef _

#define _(x) (x)

#endif

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

/*                           CompositionModes                           */

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

std::map<CompositionMode, std::string> CompositionModes()

{

    return {

        {CompositionMode::SRC_OVER, "src-over"},

        {CompositionMode::HSV_VALUE, "hsv-value"},

        {CompositionMode::MULTIPLY, "multiply"},

        {CompositionMode::SCREEN, "screen"},

        {CompositionMode::OVERLAY, "overlay"},

        {CompositionMode::HARD_LIGHT, "hard-light"},

        {CompositionMode::DARKEN, "darken"},

        {CompositionMode::LIGHTEN, "lighten"},

        {CompositionMode::COLOR_BURN, "color-burn"},

        {CompositionMode::COLOR_DODGE, "color-dodge"},

    };

}

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

/*                      CompositionModeToString()                       */

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

std::string CompositionModeToString(CompositionMode mode)

{

    const auto &modes = CompositionModes();

    const auto &iter = modes.find(mode);

    if (iter != modes.end())

    {

        return iter->second;

    }

    CPLError(CE_Failure, CPLE_IllegalArg,

             "Invalid composition mode value: %d, returning 'src-over'",

             static_cast<int>(mode));

    return "src-over";

}

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

/*                    CompositionModesIdentifiers()                     */

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

std::vector<std::string> CompositionModesIdentifiers()

{

    const auto &modes = CompositionModes();

    std::vector<std::string> identifiers;

    for (const auto &pair : modes)

    {

        identifiers.push_back(pair.second);

    }

    return identifiers;

}

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

/*                     CompositionModeFromString()                      */

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

CompositionMode CompositionModeFromString(const std::string &str)

{

    const auto &modes = CompositionModes();

    auto iter =

        std::find_if(modes.begin(), modes.end(),

                     [&str](const auto &pair) { return pair.second == str; });

    if (iter != modes.end())

    {

        return iter->first;

    }

    CPLError(CE_Failure, CPLE_IllegalArg,

             "Invalid composition identifier: %s, returning SRC_OVER",

             str.c_str());

    return CompositionMode::SRC_OVER;

}

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

/*                   MinBandCountForCompositionMode()                   */

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

//! Returns the minimum number of bands required for the given composition mode

int MinBandCountForCompositionMode(CompositionMode mode)

{

    switch (mode)

    {

        case CompositionMode::HSV_VALUE:

            return 3;

        case CompositionMode::SRC_OVER:

        case CompositionMode::MULTIPLY:

        case CompositionMode::SCREEN:

        case CompositionMode::OVERLAY:

        case CompositionMode::HARD_LIGHT:

        case CompositionMode::DARKEN:

        case CompositionMode::LIGHTEN:

        case CompositionMode::COLOR_BURN:

        case CompositionMode::COLOR_DODGE:

            return 1;

    }

    // unreachable...

    return 1;

}

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

/*                   MaxBandCountForCompositionMode()                   */

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

/**

 *  Returns the maximum number of bands allowed for the given composition mode

 */

int MaxBandCountForCompositionMode(CompositionMode mode)

{

    switch (mode)

    {

        case CompositionMode::SRC_OVER:

        case CompositionMode::HSV_VALUE:

        case CompositionMode::MULTIPLY:

        case CompositionMode::SCREEN:

        case CompositionMode::OVERLAY:

        case CompositionMode::HARD_LIGHT:

        case CompositionMode::DARKEN:

        case CompositionMode::LIGHTEN:

        case CompositionMode::COLOR_BURN:

        case CompositionMode::COLOR_DODGE:

            return 4;

    }

    // unreachable...

    return 4;

}

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

/*              BandCountIsCompatibleWithCompositionMode()              */

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

//! Checks whether the number of bands is compatible with the given composition mode

bool BandCountIsCompatibleWithCompositionMode(int bandCount,

                                              CompositionMode mode)

{

    const int minBands = MinBandCountForCompositionMode(mode);

    const int maxBands = MaxBandCountForCompositionMode(mode);

    return minBands <= bandCount && bandCount <= maxBands;

}

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

/*                            MulScale255()                             */

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

/** Multiply 2 bytes considering them as ratios with 255 = 100%, and return their product unscaled to [0, 255], by ceiling */

inline GByte MulScale255(GByte a, GByte b)

{

    return static_cast<GByte>((a * b + 255) / 256);

}

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

/*                         ProcessAlphaChannels                         */

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

static inline void ProcessAlphaChannels(size_t i,

                                        const GByte *CPL_RESTRICT pabyA,

                                        const GByte *CPL_RESTRICT pabyOverlayA,

                                        int nOpacity, bool bSwappedOpacity,

                                        GByte &outA, GByte &outOverlaA,

                                        GByte &outFinalAlpha)

{

    // Apply opacity depending on whether overlay and base were swapped

    const GByte byOpacity = static_cast<GByte>(nOpacity);

    if (!bSwappedOpacity)

    {

        outOverlaA =

            pabyOverlayA ? MulScale255(pabyOverlayA[i], byOpacity) : byOpacity;

        outA = pabyA ? pabyA[i] : 255;

    }

    else

    {

        outOverlaA = pabyOverlayA ? pabyOverlayA[i] : 255;

        if (outOverlaA != 255)

        {

            outOverlaA = pabyOverlayA ? pabyOverlayA[i] : 255;

        }

        outA = pabyA ? MulScale255(pabyA[i], byOpacity) : byOpacity;

    }

    // Da'  = Sa + Da - Sa.Da

    outFinalAlpha =

        static_cast<GByte>(outOverlaA + outA - MulScale255(outOverlaA, outA));

}

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

/*                            DivScale255()                             */

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

/** Divide 2 bytes considering them as ratios with 255 = 100%, and return their quotient unscaled to [0, 255], by flooring

 *  \warning Caution: this function does not check that the result actually fits on a byte, and just casts the computed value to byte.

 */

inline GByte DivScale255(GByte a, GByte b)

{

    if (a == 0)

    {

        return 0;

    }

    else if (b == 0)

    {

        return 255;

    }

    else

    {

        const int nRes = (a * 255) / b;

#ifdef DEBUG

        CPLAssert(nRes <= 255);

#endif

        return static_cast<GByte>(nRes);

    }

}

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

/*                         PremultiplyChannels                          */

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

//! Premultiply RGB channels by alpha (A)

static inline void PremultiplyChannels(size_t i, const GByte *pabyR,

                                       const GByte *pabyG, const GByte *pabyB,

                                       GByte &outR, GByte &outG, GByte &outB,

                                       const GByte &A)

{

    if (A == 255)

    {

        outR = pabyR ? pabyR[i] : 255;

        outG = pabyG ? pabyG[i] : outR;  // in case only R is present

        outB = pabyB ? pabyB[i] : outR;  // in case only R is present

    }

    else

    {

        outR = pabyR ? MulScale255(pabyR[i], A) : A;

        outG = pabyG ? MulScale255(pabyG[i], A)

                     : outR;  // in case only R is present

        outB = pabyB ? MulScale255(pabyB[i], A)

                     : outR;  // in case only R is present

    }

}

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

/*         GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm()         */

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

GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm(bool standaloneStep)

    : GDALRasterPipelineStepAlgorithm(

          NAME, DESCRIPTION, HELP_URL,

          ConstructorOptions()

              .SetStandaloneStep(standaloneStep)

              .SetAddDefaultArguments(false)

              .SetInputDatasetHelpMsg(_("Input raster dataset"))

              .SetInputDatasetAlias("color-input")

              .SetInputDatasetMetaVar("COLOR-INPUT")

              .SetOutputDatasetHelpMsg(_("Output raster dataset")))

{

    const auto AddOverlayDatasetArg = [this]()

    {

        auto &arg = AddArg("overlay", 0, _("Overlay dataset"),

                           &m_overlayDataset, GDAL_OF_RASTER)

                        .SetPositional()

                        .SetRequired();

        SetAutoCompleteFunctionForFilename(arg, GDAL_OF_RASTER);

    };

    if (standaloneStep)

    {

        AddRasterInputArgs(false, false);

        AddOverlayDatasetArg();

        AddProgressArg();

        AddRasterOutputArgs(false);

    }

    else

    {

        AddRasterHiddenInputDatasetArg();

        AddOverlayDatasetArg();

    }

    const std::vector<std::string> compositionModeChoices{

        CompositionModesIdentifiers()};

    AddArg("operator", 0, _("Composition operator"), &m_operatorIdentifier)

        .SetChoices(compositionModeChoices)

        .SetDefault(CompositionModeToString(CompositionMode::SRC_OVER))

        .AddAction(

            [this]()

            { m_operator = CompositionModeFromString(m_operatorIdentifier); });

    AddArg("opacity", 0,

           _("Opacity percentage to apply to the overlay dataset (0=fully "

             "transparent, 100=full use of overlay opacity)"),

           &m_opacity)

        .SetDefault(m_opacity)

        .SetMinValueIncluded(0)

        .SetMaxValueIncluded(OPACITY_INPUT_RANGE);

    AddValidationAction([this]() { return ValidateGlobal(); });

}

namespace

{

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

/*                             BlendDataset                             */

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

class BlendDataset final : public GDALDataset

{

  public:

    BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,

                 const CompositionMode eOperator, int nOpacity255Scale,

                 bool bSwappedOpacity);

    ~BlendDataset() override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override

    {

        return m_oColorDS.GetGeoTransform(gt);

    }

    const OGRSpatialReference *GetSpatialRef() const override

    {

        return m_oColorDS.GetSpatialRef();

    }

    bool AcquireSourcePixels(int nXOff, int nYOff, int nXSize, int nYSize,

                             int nBufXSize, int nBufYSize,

                             GDALRasterIOExtraArg *psExtraArg);

  protected:

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

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

                     GDALDataType eBufType, int nBandCount,

                     BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,

                     GSpacing nLineSpace, GSpacing nBandSpace,

                     GDALRasterIOExtraArg *psExtraArg) override;

  private:

    friend class BlendBand;

    GDALDataset &m_oColorDS;

    GDALDataset &m_oOverlayDS;

    const CompositionMode m_operator;

    const int m_opacity255Scale;

    std::vector<std::unique_ptr<BlendDataset>> m_apoOverviews{};

    int m_nCachedXOff = 0;

    int m_nCachedYOff = 0;

    int m_nCachedXSize = 0;

    int m_nCachedYSize = 0;

    int m_nCachedBufXSize = 0;

    int m_nCachedBufYSize = 0;

    GDALRasterIOExtraArg m_sCachedExtraArg{};

    std::vector<GByte> m_abyBuffer{};

    bool m_ioError = false;

    bool m_bSwappedOpacity = false;

};

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

/*                             rgb_to_hs()                              */

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

// rgb comes in as [r,g,b] with values in the range [0,255]. The returned

// values will be with hue and saturation in the range [0,1].

// Derived from hsv_merge.py

static void rgb_to_hs(int r, int g, int b, float *h, float *s)

{

    int minc, maxc;

    if (r <= g)

    {

        if (r <= b)

        {

            minc = r;

            maxc = std::max(g, b);

        }

        else /* b < r */

        {

            minc = b;

            maxc = g;

        }

    }

    else /* g < r */

    {

        if (g <= b)

        {

            minc = g;

            maxc = std::max(r, b);

        }

        else /* b < g */

        {

            minc = b;

            maxc = r;

        }

    }

    const int maxc_minus_minc = maxc - minc;

    if (s)

        *s = maxc_minus_minc / static_cast<float>(std::max(1, maxc));

    if (h)

    {

        const float maxc_minus_minc_times_6 =

            maxc_minus_minc == 0 ? 1.0f : 6.0f * maxc_minus_minc;

        if (maxc == b)

            *h = 4.0f / 6.0f + (r - g) / maxc_minus_minc_times_6;

        else if (maxc == g)

            *h = 2.0f / 6.0f + (b - r) / maxc_minus_minc_times_6;

        else

        {

            const float tmp = (g - b) / maxc_minus_minc_times_6;

            *h = tmp < 0.0f ? tmp + 1.0f : tmp;

        }

    }

}

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

/*                            choose_among()                            */

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

template <typename T>

static inline T choose_among(int idx, T a0, T a1, T a2, T a3, T a4, T a5)

{

    switch (idx)

    {

        case 0:

            return a0;

        case 1:

            return a1;

        case 2:

            return a2;

        case 3:

            return a3;

        case 4:

            return a4;

        default:

            break;

    }

    return a5;

}

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

/*                             hsv_to_rgb()                             */

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

// hsv comes in as [h,s,v] with hue and saturation in the range [0,1],

// but value in the range [0,255].

// Derived from hsv_merge.py

static void hsv_to_rgb(float h, float s, GByte v, GByte *r, GByte *g, GByte *b)

{

    const int i = static_cast<int>(6.0f * h);

    const float f = 6.0f * h - i;

    const GByte p = static_cast<GByte>(v * (1.0f - s) + 0.5f);

    const GByte q = static_cast<GByte>(v * (1.0f - s * f) + 0.5f);

    const GByte t = static_cast<GByte>(v * (1.0f - s * (1.0f - f)) + 0.5f);

    if (r)

        *r = choose_among(i, v, q, p, p, t, v);

    if (g)

        *g = choose_among(i, t, v, v, q, p, p);

    if (b)

        *b = choose_among(i, p, p, t, v, v, q);

}

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

/*                           XMM_RGB_to_HS()                            */

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

#ifdef HAVE_SSE2

static inline void

XMM_RGB_to_HS(const GByte *CPL_RESTRICT pInR, const GByte *CPL_RESTRICT pInG,

              const GByte *CPL_RESTRICT pInB, const XMMReg4Float &zero,

              const XMMReg4Float &one, const XMMReg4Float &six,

              const XMMReg4Float &two_over_six,

              const XMMReg4Float &four_over_six, XMMReg4Float &h,

              XMMReg4Float &s)

{

    const auto r = XMMReg4Float::Load4Val(pInR);

    const auto g = XMMReg4Float::Load4Val(pInG);

    const auto b = XMMReg4Float::Load4Val(pInB);

    const auto minc = XMMReg4Float::Min(XMMReg4Float::Min(r, g), b);

    const auto maxc = XMMReg4Float::Max(XMMReg4Float::Max(r, g), b);

    const auto max_minus_min = maxc - minc;

    s = max_minus_min / XMMReg4Float::Max(one, maxc);

    const auto inv_max_minus_min_times_6_0 =

        XMMReg4Float::Ternary(XMMReg4Float::Equals(max_minus_min, zero), one,

                              six * max_minus_min)

            .inverse();

    const auto tmp = (g - b) * inv_max_minus_min_times_6_0;

    h = XMMReg4Float::Ternary(

        XMMReg4Float::Equals(maxc, b),

        four_over_six + (r - g) * inv_max_minus_min_times_6_0,

        XMMReg4Float::Ternary(

            XMMReg4Float::Equals(maxc, g),

            two_over_six + (b - r) * inv_max_minus_min_times_6_0,

            XMMReg4Float::Ternary(XMMReg4Float::Lesser(tmp, zero), tmp + one,

                                  tmp)));

}

#endif

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

/*                          patch_value_line()                          */

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

static

#ifdef __GNUC__

    __attribute__((__noinline__))

#endif

    void patch_value_line(int nCount, const GByte *CPL_RESTRICT pInR,

                          const GByte *CPL_RESTRICT pInG,

                          const GByte *CPL_RESTRICT pInB,

                          const GByte *CPL_RESTRICT pInGray,

                          GByte *CPL_RESTRICT pOutR, GByte *CPL_RESTRICT pOutG,

                          GByte *CPL_RESTRICT pOutB)

{

    int i = 0;

#ifdef HAVE_SSE2

    const auto zero = XMMReg4Float::Zero();

    const auto one = XMMReg4Float::Set1(1.0f);

    const auto six = XMMReg4Float::Set1(6.0f);

    const auto two_over_six = XMMReg4Float::Set1(2.0f / 6.0f);

    const auto four_over_six = two_over_six + two_over_six;

    constexpr int ELTS = 8;

    for (; i + (ELTS - 1) < nCount; i += ELTS)

    {

        XMMReg4Float h0, s0;

        XMM_RGB_to_HS(pInR + i, pInG + i, pInB + i, zero, one, six,

                      two_over_six, four_over_six, h0, s0);

        XMMReg4Float h1, s1;

        XMM_RGB_to_HS(pInR + i + ELTS / 2, pInG + i + ELTS / 2,

                      pInB + i + ELTS / 2, zero, one, six, two_over_six,

                      four_over_six, h1, s1);

        XMMReg4Float v0, v1;

        XMMReg4Float::Load8Val(pInGray + i, v0, v1);

        const auto half = XMMReg4Float::Set1(0.5f);

        const auto six_h0 = six * h0;

        const auto idx0 = six_h0.truncate_to_int();

        const auto f0 = six_h0 - idx0.cast_to_float();

        const auto p0 = (v0 * (one - s0) + half).truncate_to_int();

        const auto q0 = (v0 * (one - s0 * f0) + half).truncate_to_int();

        const auto t0 = (v0 * (one - s0 * (one - f0)) + half).truncate_to_int();

        const auto six_h1 = six * h1;

        const auto idx1 = six_h1.truncate_to_int();

        const auto f1 = six_h1 - idx1.cast_to_float();

        const auto p1 = (v1 * (one - s1) + half).truncate_to_int();

        const auto q1 = (v1 * (one - s1 * f1) + half).truncate_to_int();

        const auto t1 = (v1 * (one - s1 * (one - f1)) + half).truncate_to_int();

        const auto idx = XMMReg8Byte::Pack(idx0, idx1);

        const auto v =

            XMMReg8Byte::Pack(v0.truncate_to_int(), v1.truncate_to_int());

        const auto p = XMMReg8Byte::Pack(p0, p1);

        const auto q = XMMReg8Byte::Pack(q0, q1);

        const auto t = XMMReg8Byte::Pack(t0, t1);

        const auto equalsTo0 = XMMReg8Byte::Equals(idx, XMMReg8Byte::Zero());

        const auto one8Byte = XMMReg8Byte::Set1(1);

        const auto equalsTo1 = XMMReg8Byte::Equals(idx, one8Byte);

        const auto two8Byte = one8Byte + one8Byte;

        const auto equalsTo2 = XMMReg8Byte::Equals(idx, two8Byte);

        const auto four8Byte = two8Byte + two8Byte;

        const auto equalsTo4 = XMMReg8Byte::Equals(idx, four8Byte);

        const auto equalsTo3 = XMMReg8Byte::Equals(idx, four8Byte - one8Byte);

        // clang-format off

        if (pOutR)

        {

            const auto out_r =

                XMMReg8Byte::Ternary(equalsTo0, v,

                XMMReg8Byte::Ternary(equalsTo1, q,

                XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo2, equalsTo3), p,

                XMMReg8Byte::Ternary(equalsTo4, t, v))));

            out_r.Store8Val(pOutR + i);

        }

        if (pOutG)

        {

            const auto out_g =

                XMMReg8Byte::Ternary(equalsTo0, t,

                XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo1, equalsTo2), v,

                XMMReg8Byte::Ternary(equalsTo3, q, p)));

            out_g.Store8Val(pOutG + i);

        }

        if (pOutB)

        {

            const auto out_b =

                XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo0, equalsTo1), p,

                XMMReg8Byte::Ternary(equalsTo2, t,

                XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo3, equalsTo4),

                                     v, q)));

            out_b.Store8Val(pOutB + i);

        }

        // clang-format on

    }

#endif

    for (; i < nCount; ++i)

    {

        float h, s;

        rgb_to_hs(pInR[i], pInG[i], pInB[i], &h, &s);

        hsv_to_rgb(h, s, pInGray[i], pOutR ? pOutR + i : nullptr,

                   pOutG ? pOutG + i : nullptr, pOutB ? pOutB + i : nullptr);

    }

}

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

/*                              BlendBand                               */

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

class BlendBand final : public GDALRasterBand

{

  public:

    BlendBand(BlendDataset &oBlendDataset, int nBandIn)

        : m_oBlendDataset(oBlendDataset)

    {

        nBand = nBandIn;

        nRasterXSize = oBlendDataset.GetRasterXSize();

        nRasterYSize = oBlendDataset.GetRasterYSize();

        oBlendDataset.m_oColorDS.GetRasterBand(1)->GetBlockSize(&nBlockXSize,

                                                                &nBlockYSize);

        eDataType = GDT_UInt8;

    }

    GDALColorInterp GetColorInterpretation() override

    {

        if (m_oBlendDataset.GetRasterCount() <= 2 && nBand == 1)

            return GCI_GrayIndex;

        else if (m_oBlendDataset.GetRasterCount() == 2 || nBand == 4)

            return GCI_AlphaBand;

        else

            return static_cast<GDALColorInterp>(GCI_RedBand + nBand - 1);

    }

    int GetOverviewCount() override

    {

        return static_cast<int>(m_oBlendDataset.m_apoOverviews.size());

    }

    GDALRasterBand *GetOverview(int idx) override

    {

        return idx >= 0 && idx < GetOverviewCount()

                   ? m_oBlendDataset.m_apoOverviews[idx]->GetRasterBand(nBand)

                   : nullptr;

    }

  protected:

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

    {

        int nReqXSize = 0;

        int nReqYSize = 0;

        GetActualBlockSize(nBlockXOff, nBlockYOff, &nReqXSize, &nReqYSize);

        return RasterIO(GF_Read, nBlockXOff * nBlockXSize,

                        nBlockYOff * nBlockYSize, nReqXSize, nReqYSize, pData,

                        nReqXSize, nReqYSize, GDT_UInt8, 1, nBlockXSize,

                        nullptr);

    }

    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;

  private:

    BlendDataset &m_oBlendDataset;

};

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

/*                     BlendDataset::BlendDataset()                     */

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

BlendDataset::BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,

                           const CompositionMode eOperator,

                           int nOpacity255Scale, bool bSwappedOpacity)

    : m_oColorDS(oColorDS), m_oOverlayDS(oOverlayDS), m_operator(eOperator),

      m_opacity255Scale(nOpacity255Scale), m_bSwappedOpacity(bSwappedOpacity)

{

    m_oColorDS.Reference();

    m_oOverlayDS.Reference();

    CPLAssert(oColorDS.GetRasterXSize() == oOverlayDS.GetRasterXSize());

    CPLAssert(oColorDS.GetRasterYSize() == oOverlayDS.GetRasterYSize());

    nRasterXSize = oColorDS.GetRasterXSize();

    nRasterYSize = oColorDS.GetRasterYSize();

    const int nOvrCount = oOverlayDS.GetRasterBand(1)->GetOverviewCount();

    bool bCanCreateOvr = true;

    for (int iBand = 1; iBand <= oColorDS.GetRasterCount(); ++iBand)

    {

        SetBand(iBand, std::make_unique<BlendBand>(*this, iBand));

        bCanCreateOvr =

            bCanCreateOvr &&

            (iBand > oColorDS.GetRasterCount() ||

             oColorDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount) &&

            (iBand > oOverlayDS.GetRasterCount() ||

             oOverlayDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount);

        const int nColorBandxIdx =

            iBand <= oColorDS.GetRasterCount() ? iBand : 1;

        const int nOverlayBandIdx =

            iBand <= oOverlayDS.GetRasterCount() ? iBand : 1;

        for (int iOvr = 0; iOvr < nOvrCount && bCanCreateOvr; ++iOvr)

        {

            const auto poColorOvrBand =

                oColorDS.GetRasterBand(nColorBandxIdx)->GetOverview(iOvr);

            const auto poGSOvrBand =

                oOverlayDS.GetRasterBand(nOverlayBandIdx)->GetOverview(iOvr);

            bCanCreateOvr =

                poColorOvrBand->GetDataset() != &oColorDS &&

                poColorOvrBand->GetDataset() == oColorDS.GetRasterBand(1)

                                                    ->GetOverview(iOvr)

                                                    ->GetDataset() &&

                poGSOvrBand->GetDataset() != &oOverlayDS &&

                poGSOvrBand->GetDataset() == oOverlayDS.GetRasterBand(1)

                                                 ->GetOverview(iOvr)

                                                 ->GetDataset() &&

                poColorOvrBand->GetXSize() == poGSOvrBand->GetXSize() &&

                poColorOvrBand->GetYSize() == poGSOvrBand->GetYSize();

        }

    }

    SetDescription(CPLSPrintf("Blend %s width %s", m_oColorDS.GetDescription(),

                              m_oOverlayDS.GetDescription()));

    if (nBands > 1)

    {

        SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");

    }

    if (bCanCreateOvr)

    {

        for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)

        {

            m_apoOverviews.push_back(std::make_unique<BlendDataset>(

                *(oColorDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),

                *(oOverlayDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),

                m_operator, nOpacity255Scale, bSwappedOpacity));

        }

    }

}

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

/*                    ~BlendDataset::BlendDataset()                     */

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

BlendDataset::~BlendDataset()

{

    m_oColorDS.ReleaseRef();

    m_oOverlayDS.ReleaseRef();

}

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

/*                 BlendDataset::AcquireSourcePixels()                  */

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

bool BlendDataset::AcquireSourcePixels(int nXOff, int nYOff, int nXSize,

                                       int nYSize, int nBufXSize, int nBufYSize,

                                       GDALRasterIOExtraArg *psExtraArg)

{

    if (nXOff == m_nCachedXOff && nYOff == m_nCachedYOff &&

        nXSize == m_nCachedXSize && nYSize == m_nCachedYSize &&

        nBufXSize == m_nCachedBufXSize && nBufYSize == m_nCachedBufYSize &&

        psExtraArg->eResampleAlg == m_sCachedExtraArg.eResampleAlg &&

        psExtraArg->bFloatingPointWindowValidity ==

            m_sCachedExtraArg.bFloatingPointWindowValidity &&

        (!psExtraArg->bFloatingPointWindowValidity ||

         (psExtraArg->dfXOff == m_sCachedExtraArg.dfXOff &&

          psExtraArg->dfYOff == m_sCachedExtraArg.dfYOff &&

          psExtraArg->dfXSize == m_sCachedExtraArg.dfXSize &&

          psExtraArg->dfYSize == m_sCachedExtraArg.dfYSize)))

    {

        return !m_abyBuffer.empty();

    }

    const int nColorCount = m_oColorDS.GetRasterCount();

    const int nOverlayCount = m_oOverlayDS.GetRasterCount();

    const int nCompsInBuffer = nColorCount + nOverlayCount;

    assert(nCompsInBuffer > 0);

    if (static_cast<size_t>(nBufXSize) >

        std::numeric_limits<size_t>::max() / nBufYSize / nCompsInBuffer)

    {

        CPLError(CE_Failure, CPLE_OutOfMemory,

                 "Out of memory allocating temporary buffer");

        m_abyBuffer.clear();

        m_ioError = true;

        return false;

    }

    const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;

    try

    {

        if (m_abyBuffer.size() < nPixelCount * nCompsInBuffer)

            m_abyBuffer.resize(nPixelCount * nCompsInBuffer);

    }

    catch (const std::exception &)

    {

        CPLError(CE_Failure, CPLE_OutOfMemory,

                 "Out of memory allocating temporary buffer");

        m_abyBuffer.clear();

        m_ioError = true;

        return false;

    }

    const bool bOK =

        (m_oColorDS.RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,

                             m_abyBuffer.data(), nBufXSize, nBufYSize,

                             GDT_UInt8, nColorCount, nullptr, 1, nBufXSize,

                             static_cast<GSpacing>(nPixelCount),

                             psExtraArg) == CE_None &&

         m_oOverlayDS.RasterIO(

             GF_Read, nXOff, nYOff, nXSize, nYSize,

             m_abyBuffer.data() + nPixelCount * nColorCount, nBufXSize,

             nBufYSize, GDT_UInt8, nOverlayCount, nullptr, 1, nBufXSize,

             static_cast<GSpacing>(nPixelCount), psExtraArg) == CE_None);

    if (bOK)

    {

        m_nCachedXOff = nXOff;

        m_nCachedYOff = nYOff;

        m_nCachedXSize = nXSize;

        m_nCachedYSize = nYSize;

        m_nCachedBufXSize = nBufXSize;

        m_nCachedBufYSize = nBufYSize;

        m_sCachedExtraArg = *psExtraArg;

    }

    else

    {

        m_abyBuffer.clear();

        m_ioError = true;

    }

    return bOK;

}

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

/*                             gTabInvDstA                              */

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

constexpr int SHIFT_DIV_DSTA = 8;

// Table of (255 * 256 + k/2) / k values for k in [0,255]

constexpr auto gTabInvDstA = []()

{

    std::array<uint16_t, 256> arr{};

    arr[0] = 0;

    for (int k = 1; k <= 255; ++k)

    {

        arr[k] = static_cast<uint16_t>(((255 << SHIFT_DIV_DSTA) + (k / 2)) / k);

    }

    return arr;

}();

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

/*                        BlendMultiply_Generic                         */

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

static void BlendMultiply_Generic(

    const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,

    const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,

    const GByte *CPL_RESTRICT pabyOverlayR,

    const GByte *CPL_RESTRICT pabyOverlayG,

    const GByte *CPL_RESTRICT pabyOverlayB,

    const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,

    GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,

    GByte nOpacity, int nOutputBands, bool bSwappedOpacity)

{

    // Generic formulas from Mapserver

    // Dca' = Sca.Dca + Sca.(1 - Da) + Dca.(1 - Sa)

    // Da'  = Sa + Da - Sa.Da

    // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)

    // TODO: optimize mathematically to avoid redundant computations

    GSpacing nDstOffset = 0;

    for (; i < N; ++i)

    {

        GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;

        GByte nR, nG, nB, nA;

        GByte nFinalAlpha;

        ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,

                             nA, nOverlayA, nFinalAlpha);

        PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);

        PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,

                            nOverlayR, nOverlayG, nOverlayB, nOverlayA);