/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  GDALComputedDataset and GDALComputedRasterBand
 * 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 "gdal_priv.h"
#include "vrtdataset.h"

#include <cmath>
#include <limits>

/************************************************************************/
/*                        GDALComputedDataset                           */
/************************************************************************/

class GDALComputedDataset final : public GDALDataset
{
    friend class GDALComputedRasterBand;

    const GDALComputedRasterBand::Operation m_op;
    CPLStringList m_aosOptions{};
    std::vector<std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>>
        m_bandDS{};
    std::vector<GDALRasterBand *> m_poBands{};
    VRTDataset m_oVRTDS;

    void AddSources(GDALComputedRasterBand *poBand);

    static const char *
    OperationToFunctionName(GDALComputedRasterBand::Operation op);

    GDALComputedDataset &operator=(const GDALComputedDataset &) = delete;
    GDALComputedDataset(GDALComputedDataset &&) = delete;
    GDALComputedDataset &operator=(GDALComputedDataset &&) = delete;

  public:
    GDALComputedDataset(const GDALComputedDataset &other);

    GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
                        GDALDataType eDT, int nBlockXSize, int nBlockYSize,
                        GDALComputedRasterBand::Operation op,
                        const GDALRasterBand *firstBand, double *pFirstConstant,
                        const GDALRasterBand *secondBand,
                        double *pSecondConstant);

    GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
                        GDALDataType eDT, int nBlockXSize, int nBlockYSize,
                        GDALComputedRasterBand::Operation op,
                        const std::vector<const GDALRasterBand *> &bands,
                        double constant);

    ~GDALComputedDataset() override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override
    {
        return m_oVRTDS.GetGeoTransform(gt);
    }

    const OGRSpatialReference *GetSpatialRef() const override
    {
        return m_oVRTDS.GetSpatialRef();
    }

    char **GetMetadata(const char *pszDomain) override
    {
        return m_oVRTDS.GetMetadata(pszDomain);
    }

    const char *GetMetadataItem(const char *pszName,
                                const char *pszDomain) override
    {
        return m_oVRTDS.GetMetadataItem(pszName, pszDomain);
    }

    void *GetInternalHandle(const char *pszHandleName) override
    {
        if (pszHandleName && EQUAL(pszHandleName, "VRT_DATASET"))
            return &m_oVRTDS;
        return nullptr;
    }
};

/************************************************************************/
/*                        IsComparisonOperator()                        */
/************************************************************************/

static bool IsComparisonOperator(GDALComputedRasterBand::Operation op)
{
    switch (op)
    {
        case GDALComputedRasterBand::Operation::OP_GT:
        case GDALComputedRasterBand::Operation::OP_GE:
        case GDALComputedRasterBand::Operation::OP_LT:
        case GDALComputedRasterBand::Operation::OP_LE:
        case GDALComputedRasterBand::Operation::OP_EQ:
        case GDALComputedRasterBand::Operation::OP_NE:
        case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
        case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
            return true;
        default:
            break;
    }
    return false;
}

/************************************************************************/
/*                        GDALComputedDataset()                         */
/************************************************************************/

GDALComputedDataset::GDALComputedDataset(const GDALComputedDataset &other)
    : GDALDataset(), m_op(other.m_op), m_aosOptions(other.m_aosOptions),
      m_poBands(other.m_poBands),
      m_oVRTDS(other.GetRasterXSize(), other.GetRasterYSize(),
               other.m_oVRTDS.GetBlockXSize(), other.m_oVRTDS.GetBlockYSize())
{
    nRasterXSize = other.nRasterXSize;
    nRasterYSize = other.nRasterYSize;

    auto poBand = new GDALComputedRasterBand(
        const_cast<const GDALComputedRasterBand &>(
            *cpl::down_cast<GDALComputedRasterBand *>(
                const_cast<GDALComputedDataset &>(other).GetRasterBand(1))),
        true);
    SetBand(1, poBand);

    GDALGeoTransform gt;
    if (const_cast<VRTDataset &>(other.m_oVRTDS).GetGeoTransform(gt) == CE_None)
    {
        m_oVRTDS.SetGeoTransform(gt);
    }

    if (const auto *poSRS =
            const_cast<VRTDataset &>(other.m_oVRTDS).GetSpatialRef())
    {
        m_oVRTDS.SetSpatialRef(poSRS);
    }

    m_oVRTDS.AddBand(other.m_oVRTDS.GetRasterBand(1)->GetRasterDataType(),
                     m_aosOptions.List());

    AddSources(poBand);
}

/************************************************************************/
/*                        GDALComputedDataset()                         */
/************************************************************************/

GDALComputedDataset::GDALComputedDataset(
    GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
    int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
    const GDALRasterBand *firstBand, double *pFirstConstant,
    const GDALRasterBand *secondBand, double *pSecondConstant)
    : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
{
    CPLAssert(firstBand != nullptr || secondBand != nullptr);
    if (firstBand)
        m_poBands.push_back(const_cast<GDALRasterBand *>(firstBand));
    if (secondBand)
        m_poBands.push_back(const_cast<GDALRasterBand *>(secondBand));

    nRasterXSize = nXSize;
    nRasterYSize = nYSize;

    if (auto poSrcDS = m_poBands.front()->GetDataset())
    {
        GDALGeoTransform gt;
        if (poSrcDS->GetGeoTransform(gt) == CE_None)
        {
            m_oVRTDS.SetGeoTransform(gt);
        }

        if (const auto *poSRS = poSrcDS->GetSpatialRef())
        {
            m_oVRTDS.SetSpatialRef(poSRS);
        }
    }

    if (op == GDALComputedRasterBand::Operation::OP_CAST)
    {
#ifdef DEBUG
        // Just for code coverage...
        CPL_IGNORE_RET_VAL(GDALComputedDataset::OperationToFunctionName(op));
#endif
        m_aosOptions.SetNameValue("subclass", "VRTSourcedRasterBand");
    }
    else
    {
        m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
        if (IsComparisonOperator(op))
        {
            m_aosOptions.SetNameValue("PixelFunctionType", "expression");
            if (firstBand && secondBand)
            {
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_expression",
                    CPLSPrintf(
                        "source1 %s source2",
                        GDALComputedDataset::OperationToFunctionName(op)));
            }
            else if (firstBand && pSecondConstant)
            {
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_expression",
                    CPLSPrintf("source1 %s %.17g",
                               GDALComputedDataset::OperationToFunctionName(op),
                               *pSecondConstant));
            }
            else if (pFirstConstant && secondBand)
            {
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_expression",
                    CPLSPrintf(
                        "%.17g %s source1", *pFirstConstant,
                        GDALComputedDataset::OperationToFunctionName(op)));
            }
            else
            {
                CPLAssert(false);
            }
        }
        else if (op == GDALComputedRasterBand::Operation::OP_SUBTRACT &&
                 pSecondConstant)
        {
            m_aosOptions.SetNameValue("PixelFunctionType", "sum");
            m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
                                      CPLSPrintf("%.17g", -(*pSecondConstant)));
        }
        else if (op == GDALComputedRasterBand::Operation::OP_DIVIDE)
        {
            if (pSecondConstant)
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "mul");
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_k",
                    CPLSPrintf("%.17g", 1.0 / (*pSecondConstant)));
            }
            else if (pFirstConstant)
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "inv");
                m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
                                          CPLSPrintf("%.17g", *pFirstConstant));
            }
            else
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "div");
            }
        }
        else if (op == GDALComputedRasterBand::Operation::OP_LOG)
        {
            CPLAssert(firstBand);
            CPLAssert(!secondBand);
            CPLAssert(!pFirstConstant);
            CPLAssert(!pSecondConstant);
            m_aosOptions.SetNameValue("PixelFunctionType", "expression");
            m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
                                      "log(source1)");
        }
        else if (op == GDALComputedRasterBand::Operation::OP_POW)
        {
            if (firstBand && secondBand)
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "expression");
                m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
                                          "source1 ^ source2");
            }
            else if (firstBand && pSecondConstant)
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "pow");
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_power",
                    CPLSPrintf("%.17g", *pSecondConstant));
            }
            else if (pFirstConstant && secondBand)
            {
                m_aosOptions.SetNameValue("PixelFunctionType", "exp");
                m_aosOptions.SetNameValue("_PIXELFN_ARG_base",
                                          CPLSPrintf("%.17g", *pFirstConstant));
            }
            else
            {
                CPLAssert(false);
            }
        }
        else
        {
            m_aosOptions.SetNameValue("PixelFunctionType",
                                      OperationToFunctionName(op));
            if (pSecondConstant)
                m_aosOptions.SetNameValue(
                    "_PIXELFN_ARG_k", CPLSPrintf("%.17g", *pSecondConstant));
        }
    }
    m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
    m_oVRTDS.AddBand(eDT, m_aosOptions.List());

    SetBand(1, poBand);

    AddSources(poBand);
}

/************************************************************************/
/*                        GDALComputedDataset()                         */
/************************************************************************/

GDALComputedDataset::GDALComputedDataset(
    GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
    int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
    const std::vector<const GDALRasterBand *> &bands, double constant)
    : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
{
    for (const GDALRasterBand *poIterBand : bands)
        m_poBands.push_back(const_cast<GDALRasterBand *>(poIterBand));

    nRasterXSize = nXSize;
    nRasterYSize = nYSize;

    if (auto poSrcDS = m_poBands.front()->GetDataset())
    {
        GDALGeoTransform gt;
        if (poSrcDS->GetGeoTransform(gt) == CE_None)
        {
            m_oVRTDS.SetGeoTransform(gt);
        }

        if (const auto *poSRS = poSrcDS->GetSpatialRef())
        {
            m_oVRTDS.SetSpatialRef(poSRS);
        }
    }

    m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
    if (op == GDALComputedRasterBand::Operation::OP_TERNARY)
    {
        m_aosOptions.SetNameValue("PixelFunctionType", "expression");
        m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
                                  "source1 ? source2 : source3");
    }
    else
    {
        m_aosOptions.SetNameValue("PixelFunctionType",
                                  OperationToFunctionName(op));
        if (!std::isnan(constant))
        {
            m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
                                      CPLSPrintf("%.17g", constant));
        }
        m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
    }
    m_oVRTDS.AddBand(eDT, m_aosOptions.List());

    SetBand(1, poBand);

    AddSources(poBand);
}

/************************************************************************/
/*                       ~GDALComputedDataset()                         */
/************************************************************************/

GDALComputedDataset::~GDALComputedDataset() = default;

/************************************************************************/
/*                       HaveAllBandsSameNoDataValue()                  */
/************************************************************************/

static bool HaveAllBandsSameNoDataValue(GDALRasterBand **apoBands,
                                        size_t nBands, bool &hasAtLeastOneNDV,
                                        double &singleNDV)
{
    hasAtLeastOneNDV = false;
    singleNDV = 0;

    int bFirstBandHasNoData = false;
    for (size_t i = 0; i < nBands; ++i)
    {
        int bHasNoData = false;
        const double dfNoData = apoBands[i]->GetNoDataValue(&bHasNoData);
        if (bHasNoData)
            hasAtLeastOneNDV = true;
        if (i == 0)
        {
            bFirstBandHasNoData = bHasNoData;
            singleNDV = dfNoData;
        }
        else if (bHasNoData != bFirstBandHasNoData)
        {
            return false;
        }
        else if (bFirstBandHasNoData &&
                 !((std::isnan(singleNDV) && std::isnan(dfNoData)) ||
                   (singleNDV == dfNoData)))
        {
            return false;
        }
    }
    return true;
}

/************************************************************************/
/*                  GDALComputedDataset::AddSources()                   */
/************************************************************************/

void GDALComputedDataset::AddSources(GDALComputedRasterBand *poBand)
{
    auto poSourcedRasterBand =
        cpl::down_cast<VRTSourcedRasterBand *>(m_oVRTDS.GetRasterBand(1));

    bool hasAtLeastOneNDV = false;
    double singleNDV = 0;
    const bool bSameNDV = HaveAllBandsSameNoDataValue(
        m_poBands.data(), m_poBands.size(), hasAtLeastOneNDV, singleNDV);

    // For inputs that are instances of GDALComputedDataset, clone them
    // to make sure we do not depend on temporary instances,
    // such as "a + b + c", which is evaluated as "(a + b) + c", and the
    // temporary band/dataset corresponding to a + b will go out of scope
    // quickly.
    for (GDALRasterBand *&band : m_poBands)
    {
        auto poDS = band->GetDataset();
        if (auto poComputedDS = dynamic_cast<GDALComputedDataset *>(poDS))
        {
            auto poComputedDSNew =
                std::make_unique<GDALComputedDataset>(*poComputedDS);
            band = poComputedDSNew->GetRasterBand(1);
            m_bandDS.emplace_back(poComputedDSNew.release());
        }

        int bHasNoData = false;
        const double dfNoData = band->GetNoDataValue(&bHasNoData);
        if (bHasNoData)
        {
            poSourcedRasterBand->AddComplexSource(band, -1, -1, -1, -1, -1, -1,
                                                  -1, -1, 0, 1, dfNoData);
        }
        else
        {
            poSourcedRasterBand->AddSimpleSource(band);
        }
        poSourcedRasterBand->m_papoSources.back()->SetName(CPLSPrintf(
            "source%d",
            static_cast<int>(poSourcedRasterBand->m_papoSources.size())));
    }
    if (hasAtLeastOneNDV)
    {
        poBand->m_bHasNoData = true;
        if (bSameNDV)
        {
            poBand->m_dfNoDataValue = singleNDV;
        }
        else
        {
            poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
        }
        poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue);
    }
}

/************************************************************************/
/*                       OperationToFunctionName()                      */
/************************************************************************/

/* static */ const char *GDALComputedDataset::OperationToFunctionName(
    GDALComputedRasterBand::Operation op)
{
    const char *ret = "";
    switch (op)
    {
        case GDALComputedRasterBand::Operation::OP_ADD:
            ret = "sum";
            break;
        case GDALComputedRasterBand::Operation::OP_SUBTRACT:
            ret = "diff";
            break;
        case GDALComputedRasterBand::Operation::OP_MULTIPLY:
            ret = "mul";
            break;
        case GDALComputedRasterBand::Operation::OP_DIVIDE:
            ret = "div";
            break;
        case GDALComputedRasterBand::Operation::OP_MIN:
            ret = "min";
            break;
        case GDALComputedRasterBand::Operation::OP_MAX:
            ret = "max";
            break;
        case GDALComputedRasterBand::Operation::OP_MEAN:
            ret = "mean";
            break;
        case GDALComputedRasterBand::Operation::OP_GT:
            ret = ">";
            break;
        case GDALComputedRasterBand::Operation::OP_GE:
            ret = ">=";
            break;
        case GDALComputedRasterBand::Operation::OP_LT:
            ret = "<";
            break;
        case GDALComputedRasterBand::Operation::OP_LE:
            ret = "<=";
            break;
        case GDALComputedRasterBand::Operation::OP_EQ:
            ret = "==";
            break;
        case GDALComputedRasterBand::Operation::OP_NE:
            ret = "!=";
            break;
        case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
            ret = "&&";
            break;
        case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
            ret = "||";
            break;
        case GDALComputedRasterBand::Operation::OP_CAST:
        case GDALComputedRasterBand::Operation::OP_TERNARY:
            break;
        case GDALComputedRasterBand::Operation::OP_ABS:
            ret = "mod";
            break;
        case GDALComputedRasterBand::Operation::OP_SQRT:
            ret = "sqrt";
            break;
        case GDALComputedRasterBand::Operation::OP_LOG:
            ret = "log";
            break;
        case GDALComputedRasterBand::Operation::OP_LOG10:
            ret = "log10";
            break;
        case GDALComputedRasterBand::Operation::OP_POW:
            ret = "pow";
            break;
    }
    return ret;
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(
    const GDALComputedRasterBand &other, bool)
    : GDALRasterBand()
{
    nRasterXSize = other.nRasterXSize;
    nRasterYSize = other.nRasterYSize;
    eDataType = other.eDataType;
    nBlockXSize = other.nBlockXSize;
    nBlockYSize = other.nBlockYSize;
}

//! @cond Doxygen_Suppress

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(
    Operation op, const std::vector<const GDALRasterBand *> &bands,
    double constant)
{
    CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN ||
              op == Operation::OP_MAX || op == Operation::OP_MEAN ||
              op == Operation::OP_TERNARY);

    CPLAssert(!bands.empty());
    nRasterXSize = bands[0]->GetXSize();
    nRasterYSize = bands[0]->GetYSize();
    eDataType = bands[0]->GetRasterDataType();
    for (size_t i = 1; i < bands.size(); ++i)
    {
        eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType());
    }

    bool hasAtLeastOneNDV = false;
    double singleNDV = 0;
    const bool bSameNDV =
        HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()),
                                    bands.size(), hasAtLeastOneNDV, singleNDV);

    if (!bSameNDV)
    {
        eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32;
    }
    else if (op == Operation::OP_TERNARY)
    {
        CPLAssert(bands.size() == 3);
        eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(),
                                      bands[2]->GetRasterDataType());
    }
    else if (!std::isnan(constant) && eDataType != GDT_Float64)
    {
        if (op == Operation::OP_MIN || op == Operation::OP_MAX)
        {
            eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false);
        }
        else
        {
            eDataType =
                (static_cast<double>(static_cast<float>(constant)) == constant)
                    ? GDT_Float32
                    : GDT_Float64;
        }
    }
    bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, bands, constant);
    m_poOwningDS.reset(l_poDS.release());
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
                                               const GDALRasterBand &firstBand,
                                               const GDALRasterBand &secondBand)
{
    nRasterXSize = firstBand.GetXSize();
    nRasterYSize = firstBand.GetYSize();

    bool hasAtLeastOneNDV = false;
    double singleNDV = 0;
    GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand),
                                  const_cast<GDALRasterBand *>(&secondBand)};
    const bool bSameNDV =
        HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV);

    const auto firstDT = firstBand.GetRasterDataType();
    const auto secondDT = secondBand.GetRasterDataType();
    if (!bSameNDV)
        eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64)
                        ? GDT_Float64
                        : GDT_Float32;
    else if (IsComparisonOperator(op))
        eDataType = GDT_Byte;
    else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
             secondDT == GDT_Byte)
        eDataType = GDT_UInt16;
    else if (firstDT == GDT_Float32 && secondDT == GDT_Float32)
        eDataType = GDT_Float32;
    else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) &&
             firstDT == secondDT)
        eDataType = firstDT;
    else
        eDataType = GDT_Float64;
    firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, &firstBand, nullptr, &secondBand, nullptr);
    m_poOwningDS.reset(l_poDS.release());
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant,
                                               const GDALRasterBand &band)
{
    CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op) ||
              op == Operation::OP_POW);

    nRasterXSize = band.GetXSize();
    nRasterYSize = band.GetYSize();
    const auto firstDT = band.GetRasterDataType();
    if (IsComparisonOperator(op))
        eDataType = GDT_Byte;
    else if (firstDT == GDT_Float32 &&
             static_cast<double>(static_cast<float>(constant)) == constant)
        eDataType = GDT_Float32;
    else
        eDataType = GDT_Float64;
    band.GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, nullptr, &constant, &band, nullptr);
    m_poOwningDS.reset(l_poDS.release());
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
                                               const GDALRasterBand &band,
                                               double constant)
{
    nRasterXSize = band.GetXSize();
    nRasterYSize = band.GetYSize();
    const auto firstDT = band.GetRasterDataType();
    if (IsComparisonOperator(op))
        eDataType = GDT_Byte;
    else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
             constant >= -128 && constant <= 127 &&
             std::floor(constant) == constant)
        eDataType = GDT_Byte;
    else if (firstDT == GDT_Float32 &&
             static_cast<double>(static_cast<float>(constant)) == constant)
        eDataType = GDT_Float32;
    else
        eDataType = GDT_Float64;
    band.GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, &band, nullptr, nullptr, &constant);
    m_poOwningDS.reset(l_poDS.release());
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
                                               const GDALRasterBand &band)
{
    CPLAssert(op == Operation::OP_ABS || op == Operation::OP_SQRT ||
              op == Operation::OP_LOG || op == Operation::OP_LOG10);
    nRasterXSize = band.GetXSize();
    nRasterYSize = band.GetYSize();
    eDataType =
        band.GetRasterDataType() == GDT_Float32 ? GDT_Float32 : GDT_Float64;
    band.GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, &band, nullptr, nullptr, nullptr);
    m_poOwningDS.reset(l_poDS.release());
}

/************************************************************************/
/*                       GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
                                               const GDALRasterBand &band,
                                               GDALDataType dt)
{
    CPLAssert(op == Operation::OP_CAST);
    nRasterXSize = band.GetXSize();
    nRasterYSize = band.GetYSize();
    eDataType = dt;
    band.GetBlockSize(&nBlockXSize, &nBlockYSize);
    auto l_poDS = std::make_unique<GDALComputedDataset>(
        this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
        op, &band, nullptr, nullptr, nullptr);
    m_poOwningDS.reset(l_poDS.release());
}

//! @endcond

/************************************************************************/
/*                      ~GDALComputedRasterBand()                       */
/************************************************************************/

GDALComputedRasterBand::~GDALComputedRasterBand()
{
    if (m_poOwningDS)
        cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0;
    poDS = nullptr;
}

/************************************************************************/
/*                  GDALComputedRasterBand::GetNoDataValue()            */
/************************************************************************/

double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData)
{
    if (pbHasNoData)
        *pbHasNoData = m_bHasNoData;
    return m_dfNoDataValue;
}

/************************************************************************/
/*                    GDALComputedRasterBandRelease()                   */
/************************************************************************/

/** Release a GDALComputedRasterBandH
 *
 * @since 3.12
 */
void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand)
{
    delete GDALComputedRasterBand::FromHandle(hBand);
}

/************************************************************************/
/*                           IReadBlock()                               */
/************************************************************************/

CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
                                          void *pData)
{
    auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
    return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff,
                                                        pData);
}

/************************************************************************/
/*                           IRasterIO()                                */
/************************************************************************/

CPLErr GDALComputedRasterBand::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)
{
    auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
    return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO(
        eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
        eBufType, nPixelSpace, nLineSpace, psExtraArg);
}